theme_set(theme_classic()+
          theme(text = element_text(family = "Helvetica", size = 10)))
orangeNavyPalette = c("orange","steelblue4")
df_hlm_movement <- read.csv("../data/movement_data_HLM_-14toEnd_allparticipants.csv")
names(df_hlm_movement)[6:13] <- paste0("movement_",names(df_hlm_movement)[6:13])
df_hlm_time <- read.csv("../data/time_data_HLM_-14toEnd_allparticipants.csv")
names(df_hlm_time)[6:13] <- paste0("time_",names(df_hlm_time)[6:13])
df_hlm_mood <- read.csv("../data/mood_data_HLM_-14toEnd_allparticipants.csv")

df_hlm <- Reduce(merge, list(df_hlm_movement,df_hlm_time,df_hlm_mood))

df_demo <- read.csv("../data/demo_info.csv")%>%
  mutate(Allocation = if_else(Allocation == "B","Lithium", "Control"))
  # filter(ID %in% unique(df_hlm$ID))

participant_num <- length(unique(df_hlm$ID))

df_hlm$Allocation[df_hlm$Allocation == "A"] <- "control"
df_hlm$Allocation[df_hlm$Allocation == "B"] <- "lithium"

# mood parameters from 0~4 to 1~5
df_hlm$pos_mean <- df_hlm$pos_mean+5
df_hlm$neg_mean <- df_hlm$neg_mean+5

# mood parameters need log transformation
df_hlm$posLog <- log(df_hlm$pos_mean+1)
df_hlm$negLog <- log(df_hlm$neg_mean+1)

#compute relative amplitude
df_hlm <- df_hlm%>%
  rowwise()%>%
  mutate(rel_amplitude  = (movement_M10_mean - movement_L5_mean)/(movement_M10_mean + movement_L5_mean))

#add season into df
df_hlm <- merge(df_hlm,df_demo[,c("ID","Season","Diagnosis")], by = "ID")

demo descriptive stats

#reorganize ethnicity
df_demo <- df_demo%>%
  mutate(Ethnicity = if_else(Ethnicity == "", "not reported", if_else(Ethnicity == "other brazilian" | Ethnicity == "Other croatian british","other",Ethnicity)))

df_control <- df_demo%>%filter(Allocation == "Control")
df_lithium <- df_demo%>%filter(Allocation == "Lithium")

table(df_demo[,c("Gender","Allocation")])
##       Allocation
## Gender Control Lithium
##      F       9      11
##      M       7       8
chisq.test(table(df_demo[,c("Gender","Allocation")]))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(df_demo[, c("Gender", "Allocation")])
## X-squared = 1.1504e-31, df = 1, p-value = 1
#ethnicity
table(df_demo[,c("Ethnicity","Allocation")])
##                 Allocation
## Ethnicity        Control Lithium
##   Asian                1       0
##   Black american       0       1
##   Hispanic             0       1
##   mixed                1       1
##   not reported         3       1
##   other                1       1
##   wb                   9      13
##   wo                   1       1
chisq.test(table(df_demo[,c("Ethnicity","Allocation")]))
## Warning in chisq.test(table(df_demo[, c("Ethnicity", "Allocation")])):
## Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(df_demo[, c("Ethnicity", "Allocation")])
## X-squared = 4.5032, df = 7, p-value = 0.7203
#age
df_demo%>%
  group_by(Allocation)%>%
  summarise(mean_Age = round(mean(Age),2),
            sd_Age = round(sd(Age),2))
## # A tibble: 2 × 3
##   Allocation mean_Age sd_Age
##   <chr>         <dbl>  <dbl>
## 1 Control        34.7  13.9 
## 2 Lithium        28.4   9.86
t.test(df_control$Age,df_lithium$Age)
## 
##  Welch Two Sample t-test
## 
## data:  df_control$Age and df_lithium$Age
## t = 1.5269, df = 26.485, p-value = 0.1386
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.180022 14.818180
## sample estimates:
## mean of x mean of y 
##  34.68750  28.36842
#BMI
df_demo%>%
  group_by(Allocation)%>%
  summarise(mean_BMI = round(mean(BMI),2),
            sd_BMI = round(sd(BMI),2))
## # A tibble: 2 × 3
##   Allocation mean_BMI sd_BMI
##   <chr>         <dbl>  <dbl>
## 1 Control        26.5   4.95
## 2 Lithium        25.4   6.28
t.test(df_control$BMI,df_lithium$BMI)
## 
##  Welch Two Sample t-test
## 
## data:  df_control$BMI and df_lithium$BMI
## t = 0.57297, df = 32.884, p-value = 0.5706
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.777079  4.954053
## sample estimates:
## mean of x mean of y 
##  26.49375  25.40526
#mood
df_mood <- df_hlm%>%
  filter(day >= -14 & day <= 28)%>%
  group_by(ID,Allocation)%>%
  summarise(pos_mean = round(mean(pos_mean,na.rm = TRUE),2),
            neg_mean = round(mean(neg_mean,na.rm = TRUE),2))
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
df_mood%>%
  group_by(Allocation)%>%
  summarise(mean_pos = round(mean(pos_mean),2),
            sd_po = round(sd(pos_mean),2),
            mean_neg = round(mean(neg_mean),2),
            sd_neg = round(sd(neg_mean),2))
## # A tibble: 2 × 5
##   Allocation mean_pos sd_po mean_neg sd_neg
##   <chr>         <dbl> <dbl>    <dbl>  <dbl>
## 1 control        10.4  3.01    10.7    4.18
## 2 lithium        10.7  3.5      8.71   3.25
t.test(df_mood$pos_mean[df_mood$Allocation == "lithium"],df_mood$pos_mean[df_mood$Allocation == "control"])
## 
##  Welch Two Sample t-test
## 
## data:  df_mood$pos_mean[df_mood$Allocation == "lithium"] and df_mood$pos_mean[df_mood$Allocation == "control"]
## t = 0.24204, df = 31.708, p-value = 0.8103
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.006914  2.547967
## sample estimates:
## mean of x mean of y 
##  10.71053  10.44000
t.test(df_mood$neg_mean[df_mood$Allocation == "lithium"],df_mood$neg_mean[df_mood$Allocation == "control"])
## 
##  Welch Two Sample t-test
## 
## data:  df_mood$neg_mean[df_mood$Allocation == "lithium"] and df_mood$neg_mean[df_mood$Allocation == "control"]
## t = -1.4999, df = 25.943, p-value = 0.1457
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.6624583  0.7288442
## sample estimates:
## mean of x mean of y 
##  8.710526 10.677333
#recording days
df_demo%>%
  group_by(Allocation)%>%
  summarise(median_total.recording.days = median(total.recording.days, na.rm = TRUE))
## # A tibble: 2 × 2
##   Allocation median_total.recording.days
##   <chr>                            <dbl>
## 1 Control                             50
## 2 Lithium                             50
table(df_demo[,c("total.recording.days","Allocation")])
##                     Allocation
## total.recording.days Control Lithium
##                   46       2       0
##                   47       2       1
##                   48       0       4
##                   49       2       2
##                   50       7      10
##                   52       0       1
t.test(df_lithium$total.recording.days,df_control$total.recording.days)
## 
##  Welch Two Sample t-test
## 
## data:  df_lithium$total.recording.days and df_control$total.recording.days
## t = 1.1578, df = 20.801, p-value = 0.2601
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4940028  1.7333190
## sample estimates:
## mean of x mean of y 
##  49.38889  48.76923
#missing days
df_demo%>%
  group_by(Allocation)%>%
  summarise(median_Missing.days = median(Missing.days, na.rm = TRUE))
## # A tibble: 2 × 2
##   Allocation median_Missing.days
##   <chr>                    <dbl>
## 1 Control                   17  
## 2 Lithium                   13.5
table(df_demo[,c("Missing.days","Allocation")])
##             Allocation
## Missing.days Control Lithium
##           0        1       3
##           3        1       0
##           4        1       2
##           7        0       1
##           9        0       1
##           11       1       1
##           12       1       1
##           14       1       0
##           15       0       1
##           17       2       0
##           18       0       1
##           19       1       1
##           22       2       0
##           24       1       0
##           25       0       1
##           27       0       1
##           31       0       1
##           32       0       1
##           34       1       0
##           35       0       2
t.test(df_lithium$Missing.days,df_control$Missing.days)
## 
##  Welch Two Sample t-test
## 
## data:  df_lithium$Missing.days and df_control$Missing.days
## t = 0.11917, df = 28.888, p-value = 0.906
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.598748  8.538919
## sample estimates:
## mean of x mean of y 
##  15.77778  15.30769
#BP subtype
table(df_demo[,c("Diagnosis","Allocation")])
##          Allocation
## Diagnosis Control Lithium
##       1         4       3
##       2        11      16
##       NOS       1       0
chisq.test(table(df_demo[,c("Diagnosis","Allocation")]))
## Warning in chisq.test(table(df_demo[, c("Diagnosis", "Allocation")])):
## Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(df_demo[, c("Diagnosis", "Allocation")])
## X-squared = 1.825, df = 2, p-value = 0.4015
#season
table(df_demo[,c("Season","Allocation")])
##         Allocation
## Season   Control Lithium
##   Autumn       4       4
##   Spring       3       3
##   Summer       5       3
##   Winter       4       9
chisq.test(table(df_demo[,c("Season","Allocation")]))
## Warning in chisq.test(table(df_demo[, c("Season", "Allocation")])): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(df_demo[, c("Season", "Allocation")])
## X-squared = 2.182, df = 3, p-value = 0.5355

data inspection

data availability

participant_num
## [1] 34
data_path <- "../data/"
df_hlm_movement <- read.csv(paste0(data_path,"movement_data_HLM_-14toEnd_allparticipants.csv"))
names(df_hlm_movement)[6:13] <- paste0("movement_",names(df_hlm_movement)[6:13])

df_hlm_time <- read.csv(paste0(data_path,"time_data_HLM_-14toEnd_allparticipants.csv"))
names(df_hlm_time)[6:13] <- paste0("time_",names(df_hlm_time)[6:13])

df_hlm_mood <- read.csv(paste0(data_path,"mood_data_HLM_-14toEnd_allparticipants.csv"))

df_hlm <- Reduce(merge, list(df_hlm_movement,df_hlm_time,df_hlm_mood))

df_demo <- read.csv(paste0(data_path,"/demo_info.csv"))%>%
  mutate(Allocation = if_else(Allocation == "B","Lithium", "Control"))

participant_num <- length(unique(df_hlm$ID))

df_hlm$Allocation[df_hlm$Allocation == "A"] <- "control"
df_hlm$Allocation[df_hlm$Allocation == "B"] <- "lithium"

# mood parameters from 0~4 to 1~5
df_hlm$pos_mean <- df_hlm$pos_mean+5
df_hlm$neg_mean <- df_hlm$neg_mean+5

# mood parameters need log transformation
df_hlm$posLog <- log(df_hlm$pos_mean+1)
df_hlm$negLog <- log(df_hlm$neg_mean+1)

#compute relative amplitude
df_hlm <- df_hlm%>%
  rowwise()%>%
  mutate(rel_amplitude  = (movement_M10_mean - movement_L5_mean)/(movement_M10_mean + movement_L5_mean))

#add season into df
df_hlm <- merge(df_hlm,df_demo[,c("ID","Season","Diagnosis")], by = "ID")

data availability

group_N <- count(unique(df_hlm%>%
  select(ID,Allocation)),Allocation)

# data availability by week (rolling mean)
ymax = max(group_N$n)+4
p1 <- ggplot(df_hlm%>%
  filter(!is.na(movement_M10_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "bottom",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
p2 <- ggplot(df_hlm%>%
  filter(!is.na(pos_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p3 <- ggplot( df_hlm%>%
  filter(!is.na(movement_M10_vmu))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-4, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p4 <- ggplot(df_hlm%>%
  filter(!is.na(pos_vmu))%>%
  group_by(Allocation)%>%
  count(day), aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax -4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax -4, color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())
group_N <- count(unique(df_hlm%>%
  select(ID,Allocation)),Allocation)

# data availability by week (rolling mean)
ymax = max(group_N$n)+4
p1a <- ggplot(df_hlm%>%
  filter(!is.na(movement_M10_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "bottom",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p1b <- ggplot(df_hlm%>%
  filter(!is.na(pos_mean))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood data availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p1c <- ggplot( df_hlm%>%
  filter(!is.na(movement_M10_vmu))%>%
  group_by(Allocation)%>%
  count(day),
  aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax-6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Movement model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

p1d <- ggplot(df_hlm%>%
  filter(!is.na(pos_vmu))%>%
  group_by(Allocation)%>%
  count(day), aes(x = day, y = n, color = Allocation))+
  geom_line(size=1.5)+
  geom_segment(x = -14, xend = -14, y = ymax, yend = ymax -6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  geom_segment(x = 28, xend = 28, y = ymax, yend = ymax -6, 
               color = "#999999",size=0.7, arrow =arrow(length = unit(0.3, "cm")))+
  labs(title = "Mood model output availability",
       y = "available N")+
  ylim(0,max(group_N$n)+1)+
  scale_color_manual(values = orangeNavyPalette)+
  scale_x_continuous(breaks = seq(-14,56,by = 7))+
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())
(p1a+p1b)/(p1c+p1d)+ 
  plot_annotation(tag_levels = 'A')
Data availability; grey arrows represent the period of time considered by this study.

Data availability; grey arrows represent the period of time considered by this study.

  ggsave("../figs/S1_data_availability.png",width = 8 ,height=6,dpi = 300)
df_hlm <- df_hlm%>%filter(day >= -14 & day <= 28)

density plot of movement parameters

p1d <- ggdensity(df_hlm$movement_M10_mean)+labs(title = "M10 mean",x="")+
  scale_x_continuous(breaks = seq(25,200,by = 75))
p2d <- ggdensity(df_hlm$movement_M10_s)+labs(title = "M10 noise",x="")
p3d <- ggdensity(df_hlm$movement_M10_vmu)+labs(title = "M10 volatility",x="")
p4d <- ggdensity(df_hlm$movement_M10_mu)+labs(title = "M10 mu",x="")
p5d <- ggdensity(df_hlm$movement_L5_mean)+labs(title = "L5 mean",x="")
p6d <- ggdensity(df_hlm$movement_L5_s)+labs(title = "L5 noise",x="")
p7d <- ggdensity(df_hlm$movement_L5_vmu)+labs(title = "L5 volatility",x="")
p8d <- ggdensity(df_hlm$movement_L5_mu)+labs(title = "L5 mu",x="")+
  scale_x_continuous(breaks = seq(0.1,0.5,by = 0.1))

(p1d|p2d|p3d|p4d)/(p5d|p6d|p7d|p8d)+
  plot_annotation(title = "distribution of movement parameters")
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

 ggsave("../figs/params_distribution_movement.png",width = 8 ,height=5,dpi = 300)
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

density plot of movement onset time parameters

p1d <- ggdensity(df_hlm$time_M10_mean, linetype = "dashed")+labs(title = "M10 time mean",x="")
p2d <- ggdensity(df_hlm$time_M10_s, linetype = "dashed")+labs(title = "M10 time noise",x="")
p3d <- ggdensity(df_hlm$time_M10_vmu, linetype = "dashed")+labs(title = "M10 time volatility",x="")
p4d <- ggdensity(df_hlm$time_M10_mu, linetype = "dashed")+labs(title = "M10 time mu",x="")
p5d <- ggdensity(df_hlm$time_L5_mean, linetype = "dashed")+labs(title = "L5 time mean",x="")
p6d <- ggdensity(df_hlm$time_L5_s, linetype = "dashed")+labs(title = "L5 time noise",x="")
p7d <- ggdensity(df_hlm$time_L5_vmu, linetype = "dashed")+labs(title = "L5 time volatility",x="")
p8d <- ggdensity(df_hlm$time_L5_mu, linetype = "dashed")+labs(title = "L5 time mu",x="")+
  scale_x_continuous(breaks = seq(0,1,by = 0.3))

(p1d|p2d|p3d|p4d)/(p5d|p6d|p7d|p8d)+
  plot_annotation(title = "distribution of onset time parameters")
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

 ggsave("../figs/params_distribution_onsetTime.png",width = 8 ,height=5,dpi = 300)
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Warning: Removed 695 rows containing non-finite values (`stat_density()`).
## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

density plot of mood

p1d <- ggdensity(df_hlm$posLog)+labs(title = "log PA")
p2d <- ggdensity(df_hlm$negLog)+labs(title = "log NA")
p3d <- ggdensity(df_hlm$pos_vmu)+labs(title = "PA vmu")
p4d <- ggdensity(df_hlm$neg_vmu)+labs(title = "NA vmu")

(p1d|p2d)/(p3d|p4d)
## Warning: Removed 577 rows containing non-finite values (`stat_density()`).
## Removed 577 rows containing non-finite values (`stat_density()`).
## Warning: Removed 219 rows containing non-finite values (`stat_density()`).
## Removed 219 rows containing non-finite values (`stat_density()`).

df_hlm$phase[df_hlm$day < 1] <- "pre"
df_hlm$phase[df_hlm$day >= 1 & df_hlm$day <= 7] <- "post1"
df_hlm$phase[df_hlm$day >= 8 & df_hlm$day <= 14] <- "post2"
df_hlm$phase[df_hlm$day >= 15 & df_hlm$day <= 21] <- "post3"
df_hlm$phase[df_hlm$day >= 22 & df_hlm$day <= 28] <- "post4"

df_hlm$phase <- factor(df_hlm$phase, levels = c("pre","post1","post2","post3","post4")) #,"post5","post6"

df_hlm <- df_hlm%>%
  relocate(phase, .after= day)

correlation among variables

#correlation heatmap

find.corr <- function(df,a,b){
  #df is the dataframe
  #a(IV) and b(DV) are the column index to which the model is fitted
  iv <- names(df)[a]
  dv <- names(df)[b]
  my.formula = as.formula(paste0("scale(",iv,") ~ scale(",dv,") + (1 + scale(",dv,") |ID)"))
  my.lme <- summary(lmer(my.formula,df))
  return(c(my.lme$coefficients[2,1],my.lme$coefficients[2,5])) #beta, p
}
df_hlm_corr2 <- df_hlm%>%
  select(ID,day, Allocation,movement_M10_mean,
         time_M10_mean,posLog,negLog,
         movement_M10_vmu,
         time_M10_vmu,
         pos_vmu,neg_vmu,
         movement_M10_s,
         time_M10_s,
         pos_s,neg_s)%>%
  rename(M10 = movement_M10_mean,
         timeM10 = time_M10_mean,
         logPA = posLog,
         logNA = negLog,
         timeM10_vmu = time_M10_vmu,
         M10_vmu = movement_M10_vmu,
         PA_vmu = pos_vmu,
         NA_vmu = neg_vmu,
         M10_s = movement_M10_s,
         timeM10_s = time_M10_s,
         PA_s = pos_s,
         NA_s = neg_s)

num.row2 = ncol(df_hlm_corr2)-3

my.corr.mat2 <- matrix(nrow = num.row2, ncol = num.row2)
my.corr.pmat2 <- matrix(nrow = num.row2, ncol = num.row2)

rownames(my.corr.mat2) <- names(df_hlm_corr2)[4:(num.row2+3)]
colnames(my.corr.mat2) <- names(df_hlm_corr2)[4:(num.row2+3)]
rownames(my.corr.pmat2) <- names(df_hlm_corr2)[4:(num.row2+3)]
colnames(my.corr.pmat2) <- names(df_hlm_corr2)[4:(num.row2+3)]

for (i in 1:num.row2){
  my.corr.mat2[i,i] = 1
  my.corr.pmat2[i,i] = 0
}

for (a in 4:(num.row2+3)){
  for (b in 4:(num.row2+3)){
    if (a != b){
      temp_out <- find.corr(df_hlm_corr2,a,b)
      my.corr.mat2[a-3,b-3] = temp_out[1]
      my.corr.pmat2[a-3,b-3] = temp_out[2]
    }
  }
}
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00201368 (tol = 0.002, component 1)
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
ggcorrplot::ggcorrplot(my.corr.mat2,p.mat = my.corr.pmat2,insig = "blank",colors=c("orange", "white", "navy"),lab = TRUE, digits = 1)

round(my.corr.mat2,3)
##                M10 timeM10  logPA  logNA M10_vmu timeM10_vmu PA_vmu NA_vmu
## M10          1.000  -0.031  0.194 -0.076  -0.002       0.001 -0.066 -0.029
## timeM10     -0.045   1.000 -0.061  0.047   0.059       0.071  0.003  0.041
## logPA        0.150  -0.007  1.000 -0.172   0.025       0.080  0.064  0.003
## logNA       -0.061   0.019 -0.188  1.000   0.153       0.156 -0.064  0.037
## M10_vmu     -0.007   0.037  0.058  0.178   1.000       0.592  0.283  0.331
## timeM10_vmu  0.005   0.054  0.122  0.153   0.593       1.000  0.273  0.306
## PA_vmu      -0.055   0.015  0.035 -0.054   0.143       0.036  1.000  0.245
## NA_vmu      -0.018   0.057 -0.049  0.075   0.305       0.088  0.401  1.000
## M10_s        0.022   0.038  0.040 -0.024   0.313       0.196 -0.007 -0.008
## timeM10_s   -0.030  -0.048 -0.035  0.053  -0.144      -0.162  0.045  0.005
## PA_s        -0.028   0.049  0.044  0.105   0.130       0.294 -0.053  0.077
## NA_s         0.014   0.041  0.009  0.091   0.148       0.245  0.191  0.260
##              M10_s timeM10_s   PA_s  NA_s
## M10         -0.005    -0.024  0.001 0.005
## timeM10      0.066    -0.057  0.073 0.065
## logPA        0.072    -0.060 -0.010 0.000
## logNA        0.010     0.029  0.040 0.053
## M10_vmu      0.390    -0.199  0.114 0.248
## timeM10_vmu  0.286    -0.204  0.204 0.350
## PA_vmu      -0.033    -0.035 -0.178 0.097
## NA_vmu       0.023     0.049 -0.096 0.256
## M10_s        1.000    -0.031  0.178 0.053
## timeM10_s   -0.052     1.000 -0.030 0.002
## PA_s         0.133    -0.020  1.000 0.332
## NA_s         0.083    -0.073  0.203 1.000
round(my.corr.pmat2,3)
##               M10 timeM10 logPA logNA M10_vmu timeM10_vmu PA_vmu NA_vmu M10_s
## M10         0.000   0.443 0.007 0.167   0.960       0.987  0.122  0.423 0.903
## timeM10     0.453   0.000 0.274 0.445   0.141       0.077  0.947  0.309 0.122
## logPA       0.003   0.870 0.000 0.022   0.683       0.212  0.398  0.954 0.251
## logNA       0.114   0.579 0.005 0.000   0.007       0.019  0.255  0.330 0.814
## M10_vmu     0.857   0.269 0.390 0.004   0.000       0.000  0.014  0.002 0.000
## timeM10_vmu 0.897   0.102 0.055 0.028   0.000       0.000  0.035  0.004 0.002
## PA_vmu      0.077   0.600 0.564 0.298   0.069       0.711  0.000  0.007 0.673
## NA_vmu      0.614   0.087 0.514 0.325   0.002       0.400  0.001  0.000 0.771
## M10_s       0.502   0.229 0.487 0.577   0.000       0.005  0.924  0.898 0.000
## timeM10_s   0.411   0.256 0.540 0.408   0.065       0.098  0.638  0.942 0.462
## PA_s        0.609   0.097 0.378 0.180   0.058       0.002  0.621  0.314 0.040
## NA_s        0.727   0.137 0.828 0.041   0.040       0.003  0.034  0.003 0.296
##             timeM10_s  PA_s  NA_s
## M10             0.512 0.983 0.930
## timeM10         0.242 0.140 0.209
## logPA           0.309 0.860 0.999
## logNA           0.621 0.437 0.218
## M10_vmu         0.079 0.351 0.021
## timeM10_vmu     0.095 0.158 0.002
## PA_vmu          0.682 0.102 0.296
## NA_vmu          0.577 0.384 0.032
## M10_s           0.739 0.025 0.530
## timeM10_s       0.000 0.719 0.980
## PA_s            0.791 0.000 0.000
## NA_s            0.303 0.010 0.000
ggsave("../figs/S2_corr_plot.png",width = 6 ,height=6,dpi = 300)

pre-randomisation check MANOVA

# check M10
df1 <- df_hlm%>%
  dplyr::filter(phase == "pre")%>%
  group_by(Allocation,ID)%>%
  summarise_at(vars(movement_L5_mean:time_M10_mu), ~mean(.x, na.rm = TRUE))

m01 <- manova(cbind(movement_M10_mean,movement_M10_vmu,time_M10_mean,time_M10_vmu) ~ Allocation, data =df1)
summary.aov(m01)
##  Response movement_M10_mean :
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Allocation   1    16.9   16.90  0.0476 0.8289
## Residuals   29 10300.2  355.18               
## 
##  Response movement_M10_vmu :
##             Df Sum Sq Mean Sq F value Pr(>F)
## Allocation   1 0.4642 0.46422  2.6364 0.1153
## Residuals   29 5.1064 0.17608               
## 
##  Response time_M10_mean :
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Allocation   1   0.922  0.9217  0.1868 0.6688
## Residuals   29 143.100  4.9345               
## 
##  Response time_M10_vmu :
##             Df  Sum Sq  Mean Sq F value Pr(>F)
## Allocation   1 0.00436 0.004359  0.0484 0.8273
## Residuals   29 2.60987 0.089996               
## 
## 3 observations deleted due to missingness

None of the parameters of interest were different pre-randomization.

#linear mixed model ## mood ### raw value

model_formula_mean <- posLog ~ Allocation * phase + Age + Gender + (1|ID)

df_hlm_lmer_mood_mean <- df_hlm%>%
  select(ID,Allocation,phase,day,posLog,Age,Gender,Season)

mod_mood_mean_pos <- lmer(model_formula_mean,data = df_hlm_lmer_mood_mean)
kable(round(summary(mod_mood_mean_pos)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.439 0.165 30.858 14.787 0.000
Allocationlithium 0.075 0.102 36.666 0.736 0.466
phasepost1 -0.114 0.039 843.496 -2.930 0.003
phasepost2 -0.028 0.040 845.510 -0.713 0.476
phasepost3 0.021 0.040 845.093 0.519 0.604
phasepost4 0.033 0.040 845.056 0.824 0.410
Age 0.002 0.004 29.512 0.454 0.653
GenderM -0.277 0.095 29.882 -2.904 0.007
Allocationlithium:phasepost1 0.019 0.053 844.156 0.358 0.721
Allocationlithium:phasepost2 -0.076 0.055 844.619 -1.384 0.167
Allocationlithium:phasepost3 -0.092 0.056 844.998 -1.647 0.100
Allocationlithium:phasepost4 -0.144 0.057 845.539 -2.534 0.011
lme.dscore(mod_mood_mean_pos,df_hlm_lmer_mood_mean,type = "lme4")%>%
  as.data.frame()%>%
  kable()
t df d
Allocationlithium 0.7361081 36.66568 0.2431318
phasepost1 -2.9295452 843.49614 -0.2017384
phasepost2 -0.7132416 845.50969 -0.0490577
phasepost3 0.5190409 845.09262 0.0357091
phasepost4 0.8240219 845.05637 0.0566925
Age 0.4540456 29.51158 0.1671603
GenderM -2.9044919 29.88190 -1.0626643
Allocationlithium:phasepost1 0.3577954 844.15572 0.0246294
Allocationlithium:phasepost2 -1.3839384 844.61873 -0.0952393
Allocationlithium:phasepost3 -1.6466073 844.99783 -0.1132902
Allocationlithium:phasepost4 -2.5339154 845.53869 -0.1742831
ee1 <- as.data.frame(Effect(c("phase", "Allocation"),mod_mood_mean_pos))

stat.test <- as.data.frame(round(summary(mod_mood_mean_pos)$coefficients,3))[9:12,]%>%
  rename(p_val = `Pr(>|t|)`)

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee1[ee1$phase != "pre" & ee1$Allocation == "control",]%>%
  rename(group1 = phase))%>%
  mutate(y.position = ee1$upper[1])%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))
## Joining, by = "group1"
p8 <- ggplot(ee1, aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="mood level")+ 
  stat_pvalue_manual(stat.test, label = "p_sig",vjust = -1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p8

ggsave("../figs/hlm_mood_val.png",width = 4 ,height=3,dpi = 300)

variability

model_formula <- val ~ Allocation * var_type * phase + Age + Gender + (var_type|ID)

df_hlm_lmer_mood_variability <- df_hlm%>%
  select(ID,Allocation,day,phase,pos_vmu,pos_s,neg_s,Age,Gender,posLog)%>%
  pivot_longer(c(pos_vmu,pos_s),names_pattern =  "(.*)_(.*)",names_to = c("mood_type","var_type"),values_to = "val")%>%
  mutate(var_type = if_else(var_type == "vmu","volatility","noise"))

df_hlm_lmer_mood_variability$var_type <- factor(df_hlm_lmer_mood_variability$var_type)

df_hlm_lmer_mood_variability$var_type <- relevel(df_hlm_lmer_mood_variability$var_type, ref = "volatility")

mod_pos_vmu <- lmer(model_formula,data = df_hlm_lmer_mood_variability)
kable(round(summary(mod_pos_vmu)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -3.658 0.300 54.622 -12.198 0.000
Allocationlithium -0.078 0.321 34.460 -0.244 0.809
var_typenoise 2.688 0.268 36.165 10.012 0.000
phasepost1 -0.522 0.079 2402.420 -6.614 0.000
phasepost2 -0.747 0.080 2402.793 -9.341 0.000
phasepost3 -0.706 0.081 2403.788 -8.721 0.000
phasepost4 -0.618 0.081 2403.788 -7.635 0.000
Age -0.003 0.005 29.865 -0.597 0.555
GenderM -0.310 0.117 30.016 -2.652 0.013
Allocationlithium:var_typenoise 0.052 0.359 36.037 0.144 0.886
Allocationlithium:phasepost1 0.062 0.105 2402.447 0.586 0.558
Allocationlithium:phasepost2 0.227 0.106 2402.661 2.145 0.032
Allocationlithium:phasepost3 0.554 0.107 2403.242 5.194 0.000
Allocationlithium:phasepost4 0.742 0.107 2403.242 6.964 0.000
var_typenoise:phasepost1 0.463 0.112 2403.171 4.144 0.000
var_typenoise:phasepost2 0.611 0.113 2404.255 5.403 0.000
var_typenoise:phasepost3 0.607 0.114 2407.094 5.298 0.000
var_typenoise:phasepost4 0.475 0.114 2407.094 4.144 0.000
Allocationlithium:var_typenoise:phasepost1 -0.092 0.149 2403.255 -0.616 0.538
Allocationlithium:var_typenoise:phasepost2 -0.227 0.150 2403.887 -1.515 0.130
Allocationlithium:var_typenoise:phasepost3 -0.648 0.151 2405.567 -4.301 0.000
Allocationlithium:var_typenoise:phasepost4 -0.793 0.151 2405.567 -5.259 0.000
ee2 <- as.data.frame(Effect(c("phase", "Allocation","var_type"),mod_pos_vmu))%>%
  filter(var_type == "volatility")

stat.test_pos_vmu <- as.data.frame(round(summary(mod_pos_vmu)$coefficients,3))[11:14,]

stat.test <- stat.test_pos_vmu%>%rename(p_val = 'Pr(>|t|)')
stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee2[ee2$phase != "pre" & ee2$Allocation == "lithium",]%>%
  rename(group1 = phase), by = c("group1"))%>%
  mutate(y.position = ee2$upper[1])%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))


#plot volatility first
p9 <- ggplot(ee2%>%filter(var_type == "volatility"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="mood volatility")+
  scale_y_continuous(limits = c(-5.4,-3),breaks = seq(-5.4,-3, by = 0.4))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility"),label = "p_sig",vjust = -1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p9

ggsave("../figs/mood_variability.png",width = 4 ,height=3,dpi = 300)

movement

relative amplitude

model_formula_amp <- rel_amplitude ~ Allocation * phase + posLog + pos_vmu + Age + Gender + (1|ID)
mod_relative_amplitude <- lmer(model_formula_amp,data = df_hlm)
kable(round(summary(mod_relative_amplitude)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.770 0.035 110.921 22.316 0.000
Allocationlithium -0.018 0.016 49.303 -1.129 0.264
phasepost1 0.005 0.009 538.256 0.488 0.626
phasepost2 0.014 0.010 543.395 1.367 0.172
phasepost3 0.012 0.009 543.427 1.342 0.180
phasepost4 0.018 0.010 550.841 1.703 0.089
posLog 0.037 0.009 533.222 4.341 0.000
pos_vmu -0.003 0.003 536.463 -1.281 0.201
Age 0.000 0.001 30.355 -0.851 0.402
GenderM 0.004 0.014 32.775 0.260 0.797
Allocationlithium:phasepost1 -0.008 0.014 539.910 -0.557 0.578
Allocationlithium:phasepost2 0.000 0.014 541.911 0.013 0.990
Allocationlithium:phasepost3 -0.006 0.013 545.726 -0.411 0.681
Allocationlithium:phasepost4 -0.017 0.015 553.843 -1.128 0.260
plot_model(mod_relative_amplitude, type = "pred",terms = c("phase", "Allocation"))

raw value

model_formula_mean <- val ~ Allocation * movement_type * phase + posLog + pos_vmu  + Age + Gender + (movement_type|ID)

df_hlm_lmer_movement_mean <- df_hlm%>%
  select(ID,Allocation,phase,day,movement_L5_mean,movement_M10_mean,Age,Gender,posLog, negLog, pos_vmu,Season)%>%
  pivot_longer(c(movement_L5_mean,movement_M10_mean),names_pattern =  "movement_(.*)_mean",names_to = c("movement_type"),values_to = "val")

df_hlm_lmer_movement_mean$movement_type <- factor(df_hlm_lmer_movement_mean$movement_type)
df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "M10")

mod_movement_mean_M10 <- lmer(model_formula_mean,data = df_hlm_lmer_movement_mean)
## boundary (singular) fit: see help('isSingular')
anova(mod_movement_mean_M10)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)
## Allocation                        663     663     1   31.82   2.9724  0.094400
## movement_type                   82158   82158     1   32.16 368.0827 < 2.2e-16
## phase                             332      83     4 1091.39   0.3721  0.828647
## posLog                           2120    2120     1  757.53   9.4960  0.002134
## pos_vmu                           408     408     1  804.91   1.8274  0.176816
## Age                                22      22     1  237.17   0.0972  0.755540
## Gender                            302     302     1  266.39   1.3541  0.245596
## Allocation:movement_type          595     595     1   32.13   2.6644  0.112381
## Allocation:phase                 1511     378     4 1093.41   1.6927  0.149344
## movement_type:phase               445     111     4 1094.08   0.4988  0.736646
## Allocation:movement_type:phase   1974     494     4 1096.48   2.2111  0.065843
##                                   
## Allocation                     .  
## movement_type                  ***
## phase                             
## posLog                         ** 
## pos_vmu                           
## Age                               
## Gender                            
## Allocation:movement_type          
## Allocation:phase                  
## movement_type:phase               
## Allocation:movement_type:phase .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kable(round(summary(mod_movement_mean_M10)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 44.642 6.399 171.727 6.977 0.000
Allocationlithium -0.879 5.809 42.440 -0.151 0.880
movement_typeL5 -53.929 4.636 48.866 -11.632 0.000
phasepost1 2.363 2.643 1097.183 0.894 0.372
phasepost2 5.807 2.867 1105.181 2.026 0.043
phasepost3 5.127 2.652 1108.042 1.933 0.053
phasepost4 6.989 2.963 1115.918 2.359 0.019
posLog 4.494 1.458 757.526 3.082 0.002
pos_vmu -0.613 0.454 804.907 -1.352 0.177
Age -0.017 0.054 237.166 -0.312 0.756
GenderM 1.564 1.344 266.386 1.164 0.246
Allocationlithium:movement_typeL5 0.366 6.329 52.646 0.058 0.954
Allocationlithium:phasepost1 -8.375 3.917 1103.919 -2.138 0.033
Allocationlithium:phasepost2 -8.055 3.973 1106.937 -2.028 0.043
Allocationlithium:phasepost3 -12.124 3.849 1111.723 -3.150 0.002
Allocationlithium:phasepost4 -13.768 4.187 1119.807 -3.288 0.001
movement_typeL5:phasepost1 -1.522 3.704 1095.622 -0.411 0.681
movement_typeL5:phasepost2 -5.849 3.970 1095.939 -1.473 0.141
movement_typeL5:phasepost3 -5.839 3.670 1096.865 -1.591 0.112
movement_typeL5:phasepost4 -7.781 4.053 1089.207 -1.920 0.055
Allocationlithium:movement_typeL5:phasepost1 7.817 5.474 1097.558 1.428 0.154
Allocationlithium:movement_typeL5:phasepost2 8.191 5.525 1097.678 1.482 0.139
Allocationlithium:movement_typeL5:phasepost3 12.797 5.315 1097.084 2.408 0.016
Allocationlithium:movement_typeL5:phasepost4 14.929 5.720 1091.567 2.610 0.009
df_hlm_lmer_movement_mean$movement_type <- factor(df_hlm_lmer_movement_mean$movement_type)
df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "L5")

mod_movement_mean_L5 <- lmer(model_formula_mean,data = df_hlm_lmer_movement_mean)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00627183 (tol = 0.002, component 1)
kable(round(summary(mod_movement_mean_L5)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -9.286 5.041 757.080 -1.842 0.066
Allocationlithium -0.513 2.474 1095.870 -0.207 0.836
movement_typeM10 53.929 4.637 48.839 11.630 0.000
phasepost1 0.840 2.625 1095.144 0.320 0.749
phasepost2 -0.042 2.780 1095.401 -0.015 0.988
phasepost3 -0.712 2.561 1095.736 -0.278 0.781
phasepost4 -0.793 2.795 1095.070 -0.284 0.777
posLog 4.494 1.458 757.535 3.082 0.002
pos_vmu -0.613 0.454 804.924 -1.352 0.177
Age -0.017 0.054 237.231 -0.313 0.755
GenderM 1.565 1.344 266.457 1.165 0.245
Allocationlithium:movement_typeM10 -0.366 6.330 52.616 -0.058 0.954
Allocationlithium:phasepost1 -0.557 3.825 1095.232 -0.146 0.884
Allocationlithium:phasepost2 0.135 3.848 1095.589 0.035 0.972
Allocationlithium:phasepost3 0.673 3.691 1096.204 0.182 0.855
Allocationlithium:phasepost4 1.161 3.941 1094.894 0.295 0.768
movement_typeM10:phasepost1 1.523 3.704 1095.618 0.411 0.681
movement_typeM10:phasepost2 5.849 3.970 1095.922 1.473 0.141
movement_typeM10:phasepost3 5.839 3.670 1096.843 1.591 0.112
movement_typeM10:phasepost4 7.781 4.053 1089.167 1.920 0.055
Allocationlithium:movement_typeM10:phasepost1 -7.818 5.474 1097.546 -1.428 0.154
Allocationlithium:movement_typeM10:phasepost2 -8.191 5.525 1097.659 -1.482 0.139
Allocationlithium:movement_typeM10:phasepost3 -12.797 5.315 1097.056 -2.408 0.016
Allocationlithium:movement_typeM10:phasepost4 -14.929 5.720 1091.525 -2.610 0.009
plot_model(mod_movement_mean_M10, type = "pred", terms = c("phase", "Allocation","movement_type"))+
  scale_color_manual(values=orangeNavyPalette)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.5),
         panel.grid.minor.y = element_blank())+
  labs(title = "predicted value of movement level",y="movement level")+
  scale_y_continuous(breaks = seq(-10, 70, by = 10))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

df_hlm_lmer_movement_mean$Allocation <- factor(df_hlm_lmer_movement_mean$Allocation, levels = c("lithium","control"))
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))
## Warning: Model failed to converge with 1 negative eigenvalue: -3.9e-02
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -9.800 5.068 793.682 -1.934 0.053
Allocationcontrol 0.513 2.474 1096.082 0.207 0.836
movement_typeM10 53.563 4.308 57.573 12.434 0.000
phasepost1 0.282 2.798 1095.388 0.101 0.920
phasepost2 0.094 2.675 1095.647 0.035 0.972
phasepost3 -0.039 2.641 1095.405 -0.015 0.988
phasepost4 0.369 2.773 1094.914 0.133 0.894
posLog 4.494 1.458 757.525 3.082 0.002
pos_vmu -0.613 0.454 804.906 -1.352 0.177
Age -0.017 0.054 237.165 -0.312 0.756
GenderM 1.564 1.344 266.385 1.164 0.246
Allocationcontrol:movement_typeM10 0.366 6.329 52.644 0.058 0.954
Allocationcontrol:phasepost1 0.558 3.825 1095.233 0.146 0.884
Allocationcontrol:phasepost2 -0.135 3.848 1095.589 -0.035 0.972
Allocationcontrol:phasepost3 -0.673 3.691 1096.205 -0.182 0.855
Allocationcontrol:phasepost4 -1.161 3.941 1094.893 -0.295 0.768
movement_typeM10:phasepost1 -6.295 4.027 1095.519 -1.563 0.118
movement_typeM10:phasepost2 -2.342 3.835 1096.171 -0.611 0.541
movement_typeM10:phasepost3 -6.958 3.835 1092.424 -1.814 0.070
movement_typeM10:phasepost4 -7.148 4.039 1094.587 -1.770 0.077
Allocationcontrol:movement_typeM10:phasepost1 7.817 5.474 1097.559 1.428 0.154
Allocationcontrol:movement_typeM10:phasepost2 8.191 5.525 1097.678 1.482 0.139
Allocationcontrol:movement_typeM10:phasepost3 12.797 5.315 1097.084 2.408 0.016
Allocationcontrol:movement_typeM10:phasepost4 14.929 5.720 1091.567 2.610 0.009
df_hlm_lmer_movement_mean$Allocation <- relevel(df_hlm_lmer_movement_mean$Allocation, ref = "control")

relevel to different weeks to show group contrast

df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "M10")
df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post1")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))
## boundary (singular) fit: see help('isSingular')
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 47.004 6.440 174.854 7.298 0.000
Allocationlithium -9.254 6.089 51.066 -1.520 0.135
movement_typeL5 -55.451 4.848 58.069 -11.437 0.000
phasepre -2.363 2.643 1097.183 -0.894 0.372
phasepost2 3.444 3.024 1105.717 1.139 0.255
phasepost3 2.765 2.851 1108.841 0.970 0.332
phasepost4 4.626 3.138 1117.351 1.474 0.141
posLog 4.494 1.458 757.526 3.082 0.002
pos_vmu -0.613 0.454 804.907 -1.352 0.177
Age -0.017 0.054 237.166 -0.312 0.756
GenderM 1.564 1.344 266.386 1.164 0.246
Allocationlithium:movement_typeL5 8.183 6.800 69.303 1.203 0.233
Allocationlithium:phasepre 8.375 3.917 1103.919 2.138 0.033
Allocationlithium:phasepost2 0.320 4.386 1111.321 0.073 0.942
Allocationlithium:phasepost3 -3.749 4.259 1112.941 -0.880 0.379
Allocationlithium:phasepost4 -5.393 4.546 1118.762 -1.186 0.236
movement_typeL5:phasepre 1.522 3.704 1095.622 0.411 0.681
movement_typeL5:phasepost2 -4.326 4.204 1097.141 -1.029 0.304
movement_typeL5:phasepost3 -4.316 3.939 1097.240 -1.096 0.273
movement_typeL5:phasepost4 -6.258 4.297 1090.266 -1.456 0.146
Allocationlithium:movement_typeL5:phasepre -7.817 5.474 1097.558 -1.428 0.154
Allocationlithium:movement_typeL5:phasepost2 0.374 6.075 1099.011 0.062 0.951
Allocationlithium:movement_typeL5:phasepost3 4.980 5.877 1098.714 0.847 0.397
Allocationlithium:movement_typeL5:phasepost4 7.112 6.221 1093.065 1.143 0.253
df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post2")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))
## boundary (singular) fit: see help('isSingular')
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 50.448 6.661 197.659 7.574 0.000
Allocationlithium -8.934 6.097 51.213 -1.465 0.149
movement_typeL5 -59.777 5.041 67.253 -11.859 0.000
phasepost1 -3.444 3.024 1105.717 -1.139 0.255
phasepre -5.807 2.867 1105.181 -2.026 0.043
phasepost3 -0.680 2.887 1097.231 -0.235 0.814
phasepost4 1.182 3.218 1111.685 0.367 0.713
posLog 4.494 1.458 757.526 3.082 0.002
pos_vmu -0.613 0.454 804.907 -1.352 0.177
Age -0.017 0.054 237.166 -0.312 0.756
GenderM 1.564 1.344 266.386 1.164 0.246
Allocationlithium:movement_typeL5 8.557 6.804 69.431 1.258 0.213
Allocationlithium:phasepost1 -0.320 4.386 1111.321 -0.073 0.942
Allocationlithium:phasepre 8.055 3.973 1106.937 2.028 0.043
Allocationlithium:phasepost3 -4.069 4.083 1099.931 -0.997 0.319
Allocationlithium:phasepost4 -5.713 4.480 1115.125 -1.275 0.203
movement_typeL5:phasepost1 4.326 4.204 1097.141 1.029 0.304
movement_typeL5:phasepre 5.849 3.970 1095.939 1.473 0.141
movement_typeL5:phasepost3 0.010 4.061 1095.973 0.002 0.998
movement_typeL5:phasepost4 -1.932 4.451 1093.847 -0.434 0.664
Allocationlithium:movement_typeL5:phasepost1 -0.374 6.075 1099.011 -0.062 0.951
Allocationlithium:movement_typeL5:phasepre -8.191 5.525 1097.678 -1.482 0.139
Allocationlithium:movement_typeL5:phasepost3 4.606 5.726 1096.015 0.804 0.421
Allocationlithium:movement_typeL5:phasepost4 6.738 6.179 1094.876 1.091 0.276
df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post3")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))
## boundary (singular) fit: see help('isSingular')
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 49.769 6.579 189.266 7.565 0.000
Allocationlithium -13.003 5.980 47.452 -2.174 0.035
movement_typeL5 -59.768 4.784 54.935 -12.493 0.000
phasepost2 0.680 2.887 1097.231 0.235 0.814
phasepost1 -2.765 2.851 1108.841 -0.970 0.332
phasepre -5.127 2.652 1108.042 -1.933 0.053
phasepost4 1.862 2.978 1110.217 0.625 0.532
posLog 4.494 1.458 757.526 3.082 0.002
pos_vmu -0.613 0.454 804.907 -1.352 0.177
Age -0.017 0.054 237.166 -0.312 0.756
GenderM 1.564 1.344 266.386 1.164 0.246
Allocationlithium:movement_typeL5 13.163 6.598 61.663 1.995 0.050
Allocationlithium:phasepost2 4.069 4.083 1099.931 0.997 0.319
Allocationlithium:phasepost1 3.749 4.259 1112.941 0.880 0.379
Allocationlithium:phasepre 12.124 3.849 1111.723 3.150 0.002
Allocationlithium:phasepost4 -1.644 4.239 1111.287 -0.388 0.698
movement_typeL5:phasepost2 -0.010 4.061 1095.973 -0.002 0.998
movement_typeL5:phasepost1 4.316 3.939 1097.240 1.096 0.273
movement_typeL5:phasepre 5.839 3.670 1096.865 1.591 0.112
movement_typeL5:phasepost4 -1.942 4.139 1097.203 -0.469 0.639
Allocationlithium:movement_typeL5:phasepost2 -4.606 5.726 1096.015 -0.804 0.421
Allocationlithium:movement_typeL5:phasepost1 -4.980 5.877 1098.714 -0.847 0.397
Allocationlithium:movement_typeL5:phasepre -12.797 5.315 1097.084 -2.408 0.016
Allocationlithium:movement_typeL5:phasepost4 2.132 5.903 1099.278 0.361 0.718
df_hlm_lmer_movement_mean$phase <- relevel(df_hlm_lmer_movement_mean$phase, ref = "post4")
kable(round(summary(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean))$coefficients,3))
## boundary (singular) fit: see help('isSingular')
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 51.630 6.681 201.171 7.728 0.000
Allocationlithium -14.647 6.129 52.572 -2.390 0.020
movement_typeL5 -61.710 5.035 67.324 -12.257 0.000
phasepost3 -1.862 2.978 1110.217 -0.625 0.532
phasepost2 -1.182 3.218 1111.685 -0.367 0.713
phasepost1 -4.626 3.138 1117.351 -1.474 0.141
phasepre -6.989 2.963 1115.918 -2.359 0.019
posLog 4.494 1.458 757.526 3.082 0.002
pos_vmu -0.613 0.454 804.907 -1.352 0.177
Age -0.017 0.054 237.166 -0.312 0.756
GenderM 1.564 1.344 266.386 1.164 0.246
Allocationlithium:movement_typeL5 15.295 6.869 72.649 2.227 0.029
Allocationlithium:phasepost3 1.644 4.239 1111.287 0.388 0.698
Allocationlithium:phasepost2 5.713 4.480 1115.125 1.275 0.203
Allocationlithium:phasepost1 5.393 4.546 1118.762 1.186 0.236
Allocationlithium:phasepre 13.768 4.187 1119.807 3.288 0.001
movement_typeL5:phasepost3 1.942 4.139 1097.203 0.469 0.639
movement_typeL5:phasepost2 1.932 4.451 1093.847 0.434 0.664
movement_typeL5:phasepost1 6.258 4.297 1090.266 1.456 0.146
movement_typeL5:phasepre 7.781 4.053 1089.207 1.920 0.055
Allocationlithium:movement_typeL5:phasepost3 -2.132 5.903 1099.278 -0.361 0.718
Allocationlithium:movement_typeL5:phasepost2 -6.738 6.179 1094.876 -1.091 0.276
Allocationlithium:movement_typeL5:phasepost1 -7.112 6.221 1093.065 -1.143 0.253
Allocationlithium:movement_typeL5:phasepre -14.929 5.720 1091.567 -2.610 0.009
df_hlm_lmer_movement_mean$phase <- factor(df_hlm_lmer_movement_mean$phase,
                                          levels = c("pre","post1","post2","post3","post4"))

plot_model(lmer(model_formula_mean,data = df_hlm_lmer_movement_mean),
           type = "pred",terms = c("phase","Allocation"))
## boundary (singular) fit: see help('isSingular')

# Plot Model using another method
df_hlm_lmer_movement_mean$movement_type <- relevel(df_hlm_lmer_movement_mean$movement_type, ref = "M10")
ee1 <- as.data.frame(Effect(c("phase", "Allocation","movement_type"),mod_movement_mean_M10))

stat.test_M10 <- as.data.frame(round(summary(mod_movement_mean_M10)$coefficients,3))[13:16,]
stat.test_L5 <- as.data.frame(round(summary(mod_movement_mean_L5)$coefficients,3))[13:16,]
stat.test <- rbind(stat.test_L5,stat.test_M10)%>%rename(p_val = 'Pr(>|t|)')
stat.test[1:4,"movement_type"] <- "L5" 
stat.test[5:8,"movement_type"] <- "M10" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee1[ee1$phase != "pre" & ee1$Allocation == "control",]%>%
  rename(group1 = phase))%>%
  mutate(y.position = if_else(movement_type == "L5",ee1$upper[11],if_else(movement_type == "M10", ee1$upper[1],0)))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))
## Joining, by = c("movement_type", "group1")
p4a <- ggplot(ee1, aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_grid(~movement_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="movement level")+ 
  scale_y_continuous(limits = c(-10,80),breaks = seq(-10, 80, by = 10))+
  stat_pvalue_manual(stat.test, label = "p_sig",vjust = -2)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4a

ggsave("../figs/fig2_hlm_movement_val.png",width = 6 ,height=4,dpi = 300)

p4a_M10 <- ggplot(ee1%>%filter(movement_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime movement level", x ="")+ 
  scale_y_continuous(limits = c(40,85),breaks = seq(40, 85, by = 10))+
  stat_pvalue_manual(stat.test%>%filter(movement_type == "M10"), label = "p_sig",vjust = -4)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4a_M10

variability

model_formula <- val ~ Allocation * var_type * movement_type * phase + posLog + pos_vmu + Age + Gender + (var_type * movement_type| ID)

df_hlm_lmer_movement_variability <- df_hlm%>%
  select(ID,Allocation,day,phase,movement_L5_s,movement_L5_vmu, movement_M10_s,movement_M10_vmu,Age,Gender,posLog, negLog, pos_vmu)%>%
  pivot_longer(c(movement_L5_s,movement_L5_vmu, movement_M10_s,movement_M10_vmu),names_pattern =  "movement_(.*)_(.*)",names_to = c("movement_type","var_type"),values_to = "val")%>%
  mutate(var_type = if_else(var_type == "vmu","volatility","noise"))

df_hlm_lmer_movement_variability$var_type <- factor(df_hlm_lmer_movement_variability$var_type)
df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "noise")
df_hlm_lmer_movement_variability$movement_type <- factor(df_hlm_lmer_movement_variability$movement_type)
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")

mod_movement_var_s_M10 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_s_M10)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.658 0.167 79.075 -9.919 0.000
Allocationlithium -0.001 0.130 46.126 -0.009 0.993
var_typevolatility -2.680 0.179 39.262 -15.002 0.000
movement_typeL5 0.078 0.137 45.907 0.571 0.571
phasepost1 -0.130 0.065 3379.204 -1.988 0.047
phasepost2 -0.090 0.067 3391.826 -1.348 0.178
phasepost3 -0.008 0.067 3388.243 -0.120 0.904
phasepost4 -0.009 0.068 3387.068 -0.134 0.894
posLog 0.097 0.028 3058.970 3.416 0.001
pos_vmu 0.028 0.009 2824.919 2.962 0.003
Age 0.000 0.003 28.449 -0.152 0.881
GenderM 0.127 0.073 30.311 1.743 0.091
Allocationlithium:var_typevolatility 0.013 0.240 40.105 0.054 0.957
Allocationlithium:movement_typeL5 -0.107 0.185 47.548 -0.580 0.565
var_typevolatility:movement_typeL5 -0.111 0.206 44.125 -0.537 0.594
Allocationlithium:phasepost1 0.136 0.090 3380.514 1.513 0.130
Allocationlithium:phasepost2 0.162 0.092 3383.483 1.755 0.079
Allocationlithium:phasepost3 -0.035 0.094 3387.121 -0.374 0.708
Allocationlithium:phasepost4 -0.082 0.096 3390.937 -0.854 0.393
var_typevolatility:phasepost1 -0.712 0.092 3374.134 -7.725 0.000
var_typevolatility:phasepost2 -0.922 0.094 3383.814 -9.773 0.000
var_typevolatility:phasepost3 -1.108 0.094 3381.921 -11.776 0.000
var_typevolatility:phasepost4 -1.242 0.096 3381.777 -12.929 0.000
movement_typeL5:phasepost1 -0.043 0.092 3375.088 -0.468 0.640
movement_typeL5:phasepost2 -0.055 0.094 3387.046 -0.582 0.561
movement_typeL5:phasepost3 -0.240 0.094 3384.797 -2.550 0.011
movement_typeL5:phasepost4 -0.264 0.096 3384.655 -2.746 0.006
Allocationlithium:var_typevolatility:movement_typeL5 -0.014 0.278 45.521 -0.049 0.961
Allocationlithium:var_typevolatility:phasepost1 0.314 0.127 3377.654 2.480 0.013
Allocationlithium:var_typevolatility:phasepost2 0.164 0.130 3379.757 1.259 0.208
Allocationlithium:var_typevolatility:phasepost3 0.264 0.133 3382.035 1.993 0.046
Allocationlithium:var_typevolatility:phasepost4 0.374 0.135 3385.154 2.777 0.006
Allocationlithium:movement_typeL5:phasepost1 -0.065 0.127 3379.591 -0.510 0.610
Allocationlithium:movement_typeL5:phasepost2 -0.076 0.130 3382.147 -0.579 0.562
Allocationlithium:movement_typeL5:phasepost3 0.179 0.133 3385.068 1.348 0.178
Allocationlithium:movement_typeL5:phasepost4 0.170 0.135 3388.911 1.263 0.207
var_typevolatility:movement_typeL5:phasepost1 0.184 0.130 3375.350 1.412 0.158
var_typevolatility:movement_typeL5:phasepost2 0.107 0.133 3387.487 0.806 0.421
var_typevolatility:movement_typeL5:phasepost3 0.353 0.133 3385.124 2.651 0.008
var_typevolatility:movement_typeL5:phasepost4 0.468 0.136 3384.987 3.447 0.001
Allocationlithium:var_typevolatility:movement_typeL5:phasepost1 -0.221 0.179 3379.929 -1.235 0.217
Allocationlithium:var_typevolatility:movement_typeL5:phasepost2 0.170 0.184 3382.511 0.923 0.356
Allocationlithium:var_typevolatility:movement_typeL5:phasepost3 -0.078 0.188 3385.482 -0.416 0.678
Allocationlithium:var_typevolatility:movement_typeL5:phasepost4 0.039 0.191 3389.365 0.203 0.839
df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "volatility")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")

#same model, different reference level
mod_movement_var_vmu_M10 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_vmu_M10)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -4.338 0.216 73.002 -20.109 0.000
Allocationlithium 0.012 0.224 36.448 0.052 0.958
var_typenoise 2.680 0.179 39.265 15.002 0.000
movement_typeL5 -0.033 0.171 40.109 -0.191 0.850
phasepost1 -0.842 0.065 3375.844 -12.866 0.000
phasepost2 -1.013 0.067 3381.415 -15.079 0.000
phasepost3 -1.116 0.067 3379.165 -16.691 0.000
phasepost4 -1.251 0.068 3378.095 -18.353 0.000
posLog 0.097 0.028 3058.968 3.416 0.001
pos_vmu 0.028 0.009 2824.919 2.963 0.003
Age 0.000 0.003 28.447 -0.152 0.881
GenderM 0.127 0.073 30.309 1.743 0.091
Allocationlithium:var_typenoise -0.013 0.240 40.109 -0.054 0.957
Allocationlithium:movement_typeL5 -0.121 0.230 41.052 -0.526 0.602
var_typenoise:movement_typeL5 0.111 0.206 44.121 0.537 0.594
Allocationlithium:phasepost1 0.450 0.090 3374.121 5.019 0.000
Allocationlithium:phasepost2 0.326 0.092 3375.347 3.533 0.000
Allocationlithium:phasepost3 0.229 0.094 3377.622 2.436 0.015
Allocationlithium:phasepost4 0.293 0.096 3379.557 3.059 0.002
var_typenoise:phasepost1 0.712 0.092 3374.131 7.725 0.000
var_typenoise:phasepost2 0.922 0.094 3383.811 9.773 0.000
var_typenoise:phasepost3 1.108 0.094 3381.918 11.776 0.000
var_typenoise:phasepost4 1.242 0.096 3381.774 12.929 0.000
movement_typeL5:phasepost1 0.141 0.092 3374.016 1.530 0.126
movement_typeL5:phasepost2 0.053 0.094 3383.395 0.557 0.578
movement_typeL5:phasepost3 0.113 0.094 3381.569 1.199 0.231
movement_typeL5:phasepost4 0.204 0.096 3381.452 2.127 0.034
Allocationlithium:var_typenoise:movement_typeL5 0.014 0.278 45.516 0.049 0.961
Allocationlithium:var_typenoise:phasepost1 -0.314 0.127 3377.651 -2.480 0.013
Allocationlithium:var_typenoise:phasepost2 -0.164 0.130 3379.754 -1.259 0.208
Allocationlithium:var_typenoise:phasepost3 -0.264 0.133 3382.032 -1.993 0.046
Allocationlithium:var_typenoise:phasepost4 -0.374 0.135 3385.151 -2.777 0.006
Allocationlithium:movement_typeL5:phasepost1 -0.286 0.127 3377.448 -2.257 0.024
Allocationlithium:movement_typeL5:phasepost2 0.095 0.130 3379.453 0.725 0.468
Allocationlithium:movement_typeL5:phasepost3 0.101 0.133 3381.764 0.760 0.447
Allocationlithium:movement_typeL5:phasepost4 0.209 0.135 3384.917 1.550 0.121
var_typenoise:movement_typeL5:phasepost1 -0.184 0.130 3375.347 -1.412 0.158
var_typenoise:movement_typeL5:phasepost2 -0.107 0.133 3387.484 -0.806 0.421
var_typenoise:movement_typeL5:phasepost3 -0.353 0.133 3385.121 -2.651 0.008
var_typenoise:movement_typeL5:phasepost4 -0.468 0.136 3384.984 -3.447 0.001
Allocationlithium:var_typenoise:movement_typeL5:phasepost1 0.221 0.179 3379.926 1.235 0.217
Allocationlithium:var_typenoise:movement_typeL5:phasepost2 -0.170 0.184 3382.509 -0.923 0.356
Allocationlithium:var_typenoise:movement_typeL5:phasepost3 0.078 0.188 3385.480 0.416 0.678
Allocationlithium:var_typenoise:movement_typeL5:phasepost4 -0.039 0.191 3389.363 -0.203 0.839
df_hlm_lmer_movement_variability$Allocation <- factor(df_hlm_lmer_movement_variability$Allocation, levels = c("lithium","control"))
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -4.326 0.195 80.657 -22.234 0.000
Allocationcontrol -0.012 0.224 36.461 -0.052 0.958
var_typenoise 2.667 0.161 41.203 16.600 0.000
movement_typeL5 -0.153 0.154 42.279 -0.998 0.324
phasepost1 -0.392 0.062 3378.474 -6.361 0.000
phasepost2 -0.687 0.064 3376.060 -10.730 0.000
phasepost3 -0.887 0.066 3376.527 -13.388 0.000
phasepost4 -0.959 0.067 3380.460 -14.305 0.000
posLog 0.097 0.028 3058.989 3.416 0.001
pos_vmu 0.028 0.009 2824.958 2.963 0.003
Age 0.000 0.003 28.447 -0.152 0.880
GenderM 0.127 0.073 30.310 1.743 0.091
Allocationcontrol:var_typenoise 0.013 0.240 40.121 0.054 0.957
Allocationcontrol:movement_typeL5 0.121 0.230 41.069 0.526 0.602
var_typenoise:movement_typeL5 0.124 0.186 47.331 0.667 0.508
Allocationcontrol:phasepost1 -0.450 0.090 3374.114 -5.018 0.000
Allocationcontrol:phasepost2 -0.326 0.092 3375.340 -3.533 0.000
Allocationcontrol:phasepost3 -0.229 0.094 3377.615 -2.436 0.015
Allocationcontrol:phasepost4 -0.293 0.096 3379.551 -3.059 0.002
var_typenoise:phasepost1 0.398 0.087 3381.330 4.577 0.000
var_typenoise:phasepost2 0.758 0.090 3374.860 8.416 0.000
var_typenoise:phasepost3 0.844 0.094 3382.309 9.014 0.000
var_typenoise:phasepost4 0.868 0.095 3388.568 9.178 0.000
movement_typeL5:phasepost1 -0.145 0.087 3381.050 -1.667 0.096
movement_typeL5:phasepost2 0.147 0.090 3374.709 1.634 0.102
movement_typeL5:phasepost3 0.214 0.094 3381.972 2.283 0.022
movement_typeL5:phasepost4 0.413 0.095 3388.198 4.370 0.000
Allocationcontrol:var_typenoise:movement_typeL5 -0.014 0.278 45.534 -0.049 0.961
Allocationcontrol:var_typenoise:phasepost1 0.314 0.127 3377.644 2.480 0.013
Allocationcontrol:var_typenoise:phasepost2 0.164 0.130 3379.748 1.259 0.208
Allocationcontrol:var_typenoise:phasepost3 0.264 0.133 3382.026 1.993 0.046
Allocationcontrol:var_typenoise:phasepost4 0.374 0.135 3385.145 2.777 0.006
Allocationcontrol:movement_typeL5:phasepost1 0.286 0.127 3377.442 2.257 0.024
Allocationcontrol:movement_typeL5:phasepost2 -0.095 0.130 3379.447 -0.725 0.468
Allocationcontrol:movement_typeL5:phasepost3 -0.101 0.133 3381.757 -0.760 0.447
Allocationcontrol:movement_typeL5:phasepost4 -0.209 0.135 3384.911 -1.550 0.121
var_typenoise:movement_typeL5:phasepost1 0.037 0.123 3384.588 0.303 0.762
var_typenoise:movement_typeL5:phasepost2 -0.278 0.127 3376.255 -2.180 0.029
var_typenoise:movement_typeL5:phasepost3 -0.275 0.132 3385.856 -2.077 0.038
var_typenoise:movement_typeL5:phasepost4 -0.507 0.134 3393.269 -3.792 0.000
Allocationcontrol:var_typenoise:movement_typeL5:phasepost1 -0.221 0.179 3379.920 -1.235 0.217
Allocationcontrol:var_typenoise:movement_typeL5:phasepost2 0.170 0.184 3382.502 0.923 0.356
Allocationcontrol:var_typenoise:movement_typeL5:phasepost3 -0.078 0.188 3385.474 -0.416 0.678
Allocationcontrol:var_typenoise:movement_typeL5:phasepost4 0.039 0.191 3389.357 0.203 0.839
df_hlm_lmer_movement_variability$Allocation <- relevel(df_hlm_lmer_movement_variability$Allocation, ref = "control")

#same model, different reference level
df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "noise")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "L5")

mod_movement_var_s_L5 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_s_L5)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.580 0.163 75.855 -9.703 0.000
Allocationlithium -0.108 0.120 49.484 -0.902 0.371
var_typevolatility -2.791 0.122 49.626 -22.859 0.000
movement_typeM10 -0.078 0.137 45.901 -0.571 0.571
phasepost1 -0.173 0.065 3381.659 -2.647 0.008
phasepost2 -0.145 0.067 3397.093 -2.168 0.030
phasepost3 -0.248 0.067 3393.075 -3.713 0.000
phasepost4 -0.273 0.068 3391.963 -4.009 0.000
posLog 0.097 0.028 3058.974 3.416 0.001
pos_vmu 0.028 0.009 2824.919 2.963 0.003
Age 0.000 0.003 28.450 -0.152 0.880
GenderM 0.127 0.073 30.313 1.743 0.091
Allocationlithium:var_typevolatility -0.001 0.165 51.772 -0.004 0.997
Allocationlithium:movement_typeM10 0.107 0.185 47.542 0.580 0.565
var_typevolatility:movement_typeM10 0.111 0.206 44.117 0.537 0.594
Allocationlithium:phasepost1 0.071 0.090 3384.558 0.792 0.429
Allocationlithium:phasepost2 0.086 0.092 3388.230 0.937 0.349
Allocationlithium:phasepost3 0.144 0.094 3392.444 1.528 0.127
Allocationlithium:phasepost4 0.089 0.095 3396.827 0.928 0.353
var_typevolatility:phasepost1 -0.528 0.092 3378.822 -5.730 0.000
var_typevolatility:phasepost2 -0.815 0.094 3395.366 -8.647 0.000
var_typevolatility:phasepost3 -0.755 0.094 3392.150 -8.038 0.000
var_typevolatility:phasepost4 -0.774 0.096 3391.974 -8.067 0.000
movement_typeM10:phasepost1 0.043 0.092 3375.091 0.468 0.640
movement_typeM10:phasepost2 0.055 0.094 3387.049 0.582 0.561
movement_typeM10:phasepost3 0.240 0.094 3384.800 2.550 0.011
movement_typeM10:phasepost4 0.264 0.096 3384.657 2.746 0.006
Allocationlithium:var_typevolatility:movement_typeM10 0.014 0.278 45.512 0.049 0.961
Allocationlithium:var_typevolatility:phasepost1 0.093 0.127 3385.508 0.734 0.463
Allocationlithium:var_typevolatility:phasepost2 0.335 0.130 3389.064 2.566 0.010
Allocationlithium:var_typevolatility:phasepost3 0.187 0.133 3392.905 1.407 0.160
Allocationlithium:var_typevolatility:phasepost4 0.413 0.135 3397.488 3.069 0.002
Allocationlithium:movement_typeM10:phasepost1 0.065 0.127 3379.594 0.510 0.610
Allocationlithium:movement_typeM10:phasepost2 0.076 0.130 3382.150 0.579 0.562
Allocationlithium:movement_typeM10:phasepost3 -0.179 0.133 3385.071 -1.348 0.178
Allocationlithium:movement_typeM10:phasepost4 -0.170 0.135 3388.913 -1.263 0.207
var_typevolatility:movement_typeM10:phasepost1 -0.184 0.130 3375.353 -1.412 0.158
var_typevolatility:movement_typeM10:phasepost2 -0.107 0.133 3387.489 -0.806 0.421
var_typevolatility:movement_typeM10:phasepost3 -0.353 0.133 3385.126 -2.651 0.008
var_typevolatility:movement_typeM10:phasepost4 -0.468 0.136 3384.989 -3.447 0.001
Allocationlithium:var_typevolatility:movement_typeM10:phasepost1 0.221 0.179 3379.932 1.235 0.217
Allocationlithium:var_typevolatility:movement_typeM10:phasepost2 -0.170 0.184 3382.514 -0.923 0.356
Allocationlithium:var_typevolatility:movement_typeM10:phasepost3 0.078 0.188 3385.485 0.416 0.678
Allocationlithium:var_typevolatility:movement_typeM10:phasepost4 -0.039 0.191 3389.367 -0.203 0.839
df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "volatility")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "L5")

mod_movement_var_vmu_L5 <- lmer(model_formula,data = df_hlm_lmer_movement_variability)
kable(round(summary(mod_movement_var_vmu_L5)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -4.370 0.188 82.911 -23.216 0.000
Allocationlithium -0.109 0.175 39.862 -0.624 0.537
var_typenoise 2.791 0.122 49.639 22.861 0.000
movement_typeM10 0.033 0.171 40.112 0.191 0.850
phasepost1 -0.701 0.065 3377.827 -10.713 0.000
phasepost2 -0.960 0.067 3388.169 -14.305 0.000
phasepost3 -1.003 0.067 3384.895 -15.012 0.000
phasepost4 -1.047 0.068 3383.872 -15.365 0.000
posLog 0.097 0.028 3058.972 3.416 0.001
pos_vmu 0.028 0.009 2824.920 2.963 0.003
Age 0.000 0.003 28.449 -0.152 0.880
GenderM 0.127 0.073 30.312 1.743 0.091
Allocationlithium:var_typenoise 0.001 0.165 51.786 0.004 0.997
Allocationlithium:movement_typeM10 0.121 0.230 41.055 0.526 0.602
var_typenoise:movement_typeM10 -0.111 0.206 44.124 -0.537 0.594
Allocationlithium:phasepost1 0.164 0.090 3378.038 1.829 0.068
Allocationlithium:phasepost2 0.421 0.092 3380.176 4.559 0.000
Allocationlithium:phasepost3 0.330 0.094 3383.670 3.509 0.000
Allocationlithium:phasepost4 0.502 0.096 3387.088 5.246 0.000
var_typenoise:phasepost1 0.528 0.092 3378.819 5.730 0.000
var_typenoise:phasepost2 0.815 0.094 3395.363 8.647 0.000
var_typenoise:phasepost3 0.755 0.094 3392.147 8.038 0.000
var_typenoise:phasepost4 0.774 0.096 3391.971 8.067 0.000
movement_typeM10:phasepost1 -0.141 0.092 3374.017 -1.529 0.126
movement_typeM10:phasepost2 -0.053 0.094 3383.396 -0.557 0.578
movement_typeM10:phasepost3 -0.113 0.094 3381.571 -1.199 0.231
movement_typeM10:phasepost4 -0.204 0.096 3381.453 -2.127 0.034
Allocationlithium:var_typenoise:movement_typeM10 -0.014 0.278 45.520 -0.049 0.961
Allocationlithium:var_typenoise:phasepost1 -0.093 0.127 3385.505 -0.734 0.463
Allocationlithium:var_typenoise:phasepost2 -0.335 0.130 3389.061 -2.566 0.010
Allocationlithium:var_typenoise:phasepost3 -0.187 0.133 3392.902 -1.407 0.160
Allocationlithium:var_typenoise:phasepost4 -0.413 0.135 3397.485 -3.069 0.002
Allocationlithium:movement_typeM10:phasepost1 0.286 0.127 3377.449 2.257 0.024
Allocationlithium:movement_typeM10:phasepost2 -0.095 0.130 3379.455 -0.725 0.468
Allocationlithium:movement_typeM10:phasepost3 -0.101 0.133 3381.765 -0.760 0.447
Allocationlithium:movement_typeM10:phasepost4 -0.209 0.135 3384.919 -1.550 0.121
var_typenoise:movement_typeM10:phasepost1 0.184 0.130 3375.349 1.412 0.158
var_typenoise:movement_typeM10:phasepost2 0.107 0.133 3387.486 0.806 0.421
var_typenoise:movement_typeM10:phasepost3 0.353 0.133 3385.123 2.651 0.008
var_typenoise:movement_typeM10:phasepost4 0.468 0.136 3384.986 3.447 0.001
Allocationlithium:var_typenoise:movement_typeM10:phasepost1 -0.221 0.179 3379.928 -1.235 0.217
Allocationlithium:var_typenoise:movement_typeM10:phasepost2 0.170 0.184 3382.511 0.923 0.356
Allocationlithium:var_typenoise:movement_typeM10:phasepost3 -0.078 0.188 3385.482 -0.416 0.678
Allocationlithium:var_typenoise:movement_typeM10:phasepost4 0.039 0.191 3389.365 0.203 0.839
plot_model(mod_movement_var_s_M10, type = "pred", terms = c("phase", "Allocation","movement_type","var_type"))

Conclusion: lithium increased volatility of M10 and L5 by week 4. Conclusion: Lithium increases daytime volatility and decreases nighttime volatility of time

relevel to different weeks to show group contrast

df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")
df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post1")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -5.180 0.216 73.680 -23.959 0.000
Allocationlithium 0.462 0.225 36.779 2.052 0.047
var_typenoise 3.392 0.180 40.332 18.863 0.000
movement_typeL5 0.108 0.172 41.306 0.631 0.532
phasepre 0.842 0.065 3375.844 12.866 0.000
phasepost2 -0.171 0.068 3376.764 -2.501 0.012
phasepost3 -0.274 0.068 3376.823 -4.014 0.000
phasepost4 -0.409 0.070 3377.550 -5.858 0.000
posLog 0.097 0.028 3058.986 3.416 0.001
pos_vmu 0.028 0.009 2824.942 2.962 0.003
Age 0.000 0.003 28.451 -0.152 0.881
GenderM 0.127 0.073 30.313 1.743 0.091
Allocationlithium:var_typenoise -0.327 0.241 40.737 -1.357 0.182
Allocationlithium:movement_typeL5 -0.407 0.231 41.759 -1.764 0.085
var_typenoise:movement_typeL5 -0.073 0.208 45.936 -0.353 0.726
Allocationlithium:phasepre -0.450 0.090 3374.121 -5.019 0.000
Allocationlithium:phasepost2 -0.123 0.094 3375.672 -1.314 0.189
Allocationlithium:phasepost3 -0.220 0.096 3376.613 -2.307 0.021
Allocationlithium:phasepost4 -0.157 0.097 3376.411 -1.622 0.105
var_typenoise:phasepre -0.712 0.092 3374.131 -7.725 0.000
var_typenoise:phasepost2 0.210 0.096 3381.885 2.183 0.029
var_typenoise:phasepost3 0.396 0.096 3381.518 4.111 0.000
var_typenoise:phasepost4 0.530 0.099 3382.890 5.380 0.000
movement_typeL5:phasepre -0.141 0.092 3374.016 -1.529 0.126
movement_typeL5:phasepost2 -0.088 0.096 3381.553 -0.917 0.359
movement_typeL5:phasepost3 -0.028 0.096 3381.224 -0.292 0.770
movement_typeL5:phasepost4 0.063 0.099 3382.585 0.643 0.520
Allocationlithium:var_typenoise:movement_typeL5 0.235 0.279 46.623 0.841 0.405
Allocationlithium:var_typenoise:phasepre 0.314 0.127 3377.651 2.480 0.013
Allocationlithium:var_typenoise:phasepost2 0.150 0.133 3380.229 1.129 0.259
Allocationlithium:var_typenoise:phasepost3 0.050 0.135 3380.805 0.369 0.712
Allocationlithium:var_typenoise:phasepost4 -0.060 0.136 3380.010 -0.440 0.660
Allocationlithium:movement_typeL5:phasepre 0.286 0.127 3377.448 2.257 0.024
Allocationlithium:movement_typeL5:phasepost2 0.381 0.133 3379.933 2.866 0.004
Allocationlithium:movement_typeL5:phasepost3 0.387 0.135 3380.557 2.869 0.004
Allocationlithium:movement_typeL5:phasepost4 0.495 0.136 3379.924 3.626 0.000
var_typenoise:movement_typeL5:phasepre 0.184 0.130 3375.347 1.412 0.158
var_typenoise:movement_typeL5:phasepost2 0.077 0.136 3385.013 0.562 0.574
var_typenoise:movement_typeL5:phasepost3 -0.169 0.136 3384.586 -1.238 0.216
var_typenoise:movement_typeL5:phasepost4 -0.284 0.139 3386.394 -2.039 0.042
Allocationlithium:var_typenoise:movement_typeL5:phasepre -0.221 0.179 3379.927 -1.235 0.217
Allocationlithium:var_typenoise:movement_typeL5:phasepost2 -0.391 0.188 3383.024 -2.086 0.037
Allocationlithium:var_typenoise:movement_typeL5:phasepost3 -0.143 0.191 3383.830 -0.752 0.452
Allocationlithium:var_typenoise:movement_typeL5:phasepost4 -0.260 0.193 3383.081 -1.348 0.178
df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post2")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -5.350 0.218 75.992 -24.553 0.000
Allocationlithium 0.338 0.226 37.652 1.494 0.144
var_typenoise 3.602 0.181 41.396 19.897 0.000
movement_typeL5 0.020 0.173 42.502 0.116 0.908
phasepost1 0.171 0.068 3376.764 2.501 0.012
phasepre 1.013 0.067 3381.416 15.078 0.000
phasepost3 -0.103 0.069 3371.569 -1.504 0.133
phasepost4 -0.238 0.070 3373.140 -3.389 0.001
posLog 0.097 0.028 3058.975 3.416 0.001
pos_vmu 0.028 0.009 2824.933 2.963 0.003
Age 0.000 0.003 28.446 -0.152 0.880
GenderM 0.127 0.073 30.308 1.743 0.091
Allocationlithium:var_typenoise -0.177 0.244 42.429 -0.728 0.471
Allocationlithium:movement_typeL5 -0.026 0.233 43.661 -0.112 0.912
var_typenoise:movement_typeL5 0.003 0.210 47.683 0.015 0.988
Allocationlithium:phasepost1 0.123 0.094 3375.673 1.314 0.189
Allocationlithium:phasepre -0.326 0.092 3375.348 -3.533 0.000
Allocationlithium:phasepost3 -0.097 0.098 3373.944 -0.994 0.320
Allocationlithium:phasepost4 -0.034 0.099 3375.862 -0.338 0.735
var_typenoise:phasepost1 -0.210 0.096 3381.886 -2.183 0.029
var_typenoise:phasepre -0.922 0.094 3383.812 -9.773 0.000
var_typenoise:phasepost3 0.186 0.097 3372.000 1.912 0.056
var_typenoise:phasepost4 0.320 0.099 3374.944 3.216 0.001
movement_typeL5:phasepost1 0.088 0.096 3381.554 0.917 0.359
movement_typeL5:phasepre -0.053 0.094 3383.396 -0.557 0.578
movement_typeL5:phasepost3 0.060 0.097 3371.959 0.621 0.535
movement_typeL5:phasepost4 0.152 0.099 3374.797 1.527 0.127
Allocationlithium:var_typenoise:movement_typeL5 -0.157 0.284 49.457 -0.552 0.583
Allocationlithium:var_typenoise:phasepost1 -0.150 0.133 3380.230 -1.129 0.259
Allocationlithium:var_typenoise:phasepre 0.164 0.130 3379.755 1.259 0.208
Allocationlithium:var_typenoise:phasepost3 -0.100 0.138 3376.294 -0.727 0.467
Allocationlithium:var_typenoise:phasepost4 -0.210 0.140 3379.967 -1.502 0.133
Allocationlithium:movement_typeL5:phasepost1 -0.381 0.133 3379.934 -2.866 0.004
Allocationlithium:movement_typeL5:phasepre -0.095 0.130 3379.454 -0.725 0.468
Allocationlithium:movement_typeL5:phasepost3 0.006 0.138 3376.122 0.045 0.964
Allocationlithium:movement_typeL5:phasepost4 0.114 0.140 3379.764 0.817 0.414
var_typenoise:movement_typeL5:phasepost1 -0.077 0.136 3385.014 -0.562 0.574
var_typenoise:movement_typeL5:phasepre 0.107 0.133 3387.486 0.806 0.421
var_typenoise:movement_typeL5:phasepost3 -0.245 0.137 3372.430 -1.786 0.074
var_typenoise:movement_typeL5:phasepost4 -0.361 0.141 3376.433 -2.566 0.010
Allocationlithium:var_typenoise:movement_typeL5:phasepost1 0.391 0.188 3383.025 2.086 0.037
Allocationlithium:var_typenoise:movement_typeL5:phasepre 0.170 0.184 3382.510 0.923 0.356
Allocationlithium:var_typenoise:movement_typeL5:phasepost3 0.248 0.195 3378.150 1.274 0.203
Allocationlithium:var_typenoise:movement_typeL5:phasepost4 0.132 0.198 3382.951 0.665 0.506
df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post3")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00266921 (tol = 0.002, component 1)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -5.454 0.218 76.348 -25.001 0.000
Allocationlithium 0.241 0.227 38.148 1.062 0.295
var_typenoise 3.788 0.181 41.355 20.929 0.000
movement_typeL5 0.080 0.173 42.452 0.464 0.645
phasepost2 0.103 0.069 3371.568 1.504 0.133
phasepost1 0.274 0.068 3376.823 4.014 0.000
phasepre 1.116 0.067 3379.165 16.691 0.000
phasepost4 -0.135 0.070 3371.847 -1.926 0.054
posLog 0.097 0.028 3058.986 3.416 0.001
pos_vmu 0.028 0.009 2824.941 2.962 0.003
Age 0.000 0.003 28.450 -0.152 0.881
GenderM 0.127 0.073 30.312 1.743 0.091
Allocationlithium:var_typenoise -0.277 0.245 43.274 -1.133 0.263
Allocationlithium:movement_typeL5 -0.020 0.234 44.607 -0.085 0.933
var_typenoise:movement_typeL5 -0.242 0.210 47.611 -1.152 0.255
Allocationlithium:phasepost2 0.097 0.098 3373.943 0.994 0.320
Allocationlithium:phasepost1 0.220 0.096 3376.613 2.307 0.021
Allocationlithium:phasepre -0.229 0.094 3377.622 -2.436 0.015
Allocationlithium:phasepost4 0.063 0.100 3373.835 0.635 0.525
var_typenoise:phasepost2 -0.186 0.097 3372.000 -1.912 0.056
var_typenoise:phasepost1 -0.396 0.096 3381.518 -4.111 0.000
var_typenoise:phasepre -1.108 0.094 3381.919 -11.776 0.000
var_typenoise:phasepost4 0.134 0.099 3372.954 1.352 0.177
movement_typeL5:phasepost2 -0.060 0.097 3371.958 -0.621 0.535
movement_typeL5:phasepost1 0.028 0.096 3381.224 0.292 0.770
movement_typeL5:phasepre -0.113 0.094 3381.570 -1.199 0.231
movement_typeL5:phasepost4 0.091 0.099 3372.877 0.923 0.356
Allocationlithium:var_typenoise:movement_typeL5 0.092 0.286 50.888 0.321 0.750
Allocationlithium:var_typenoise:phasepost2 0.100 0.138 3376.293 0.727 0.467
Allocationlithium:var_typenoise:phasepost1 -0.050 0.135 3380.805 -0.369 0.712
Allocationlithium:var_typenoise:phasepre 0.264 0.133 3382.032 1.993 0.046
Allocationlithium:var_typenoise:phasepost4 -0.110 0.141 3376.995 -0.778 0.437
Allocationlithium:movement_typeL5:phasepost2 -0.006 0.138 3376.121 -0.045 0.964
Allocationlithium:movement_typeL5:phasepost1 -0.387 0.135 3380.558 -2.869 0.004
Allocationlithium:movement_typeL5:phasepre -0.101 0.133 3381.764 -0.760 0.447
Allocationlithium:movement_typeL5:phasepost4 0.108 0.141 3376.777 0.765 0.444
var_typenoise:movement_typeL5:phasepost2 0.245 0.137 3372.430 1.786 0.074
var_typenoise:movement_typeL5:phasepost1 0.169 0.136 3384.587 1.238 0.216
var_typenoise:movement_typeL5:phasepre 0.353 0.133 3385.122 2.651 0.008
var_typenoise:movement_typeL5:phasepost4 -0.115 0.140 3373.736 -0.823 0.411
Allocationlithium:var_typenoise:movement_typeL5:phasepost2 -0.248 0.195 3378.149 -1.274 0.203
Allocationlithium:var_typenoise:movement_typeL5:phasepost1 0.143 0.191 3383.831 0.752 0.452
Allocationlithium:var_typenoise:movement_typeL5:phasepre -0.078 0.188 3385.480 -0.416 0.678
Allocationlithium:var_typenoise:movement_typeL5:phasepost4 -0.117 0.200 3378.935 -0.584 0.559
df_hlm_lmer_movement_variability$phase <- relevel(df_hlm_lmer_movement_variability$phase, ref = "post4")
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_movement_variability))$coefficients,3))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00341248 (tol = 0.002, component 1)
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -5.589 0.219 76.641 -25.568 0.000
Allocationlithium 0.304 0.228 38.323 1.338 0.189
var_typenoise 3.922 0.182 42.263 21.525 0.000
movement_typeL5 0.172 0.174 43.466 0.985 0.330
phasepost3 0.135 0.070 3371.897 1.927 0.054
phasepost2 0.238 0.070 3373.189 3.389 0.001
phasepost1 0.409 0.070 3377.596 5.858 0.000
phasepre 1.251 0.068 3378.141 18.353 0.000
posLog 0.097 0.028 3059.012 3.416 0.001
pos_vmu 0.028 0.009 2824.966 2.963 0.003
Age 0.000 0.003 28.435 -0.152 0.880
GenderM 0.127 0.073 30.296 1.743 0.091
Allocationlithium:var_typenoise -0.387 0.246 43.683 -1.576 0.122
Allocationlithium:movement_typeL5 0.088 0.235 45.061 0.375 0.710
var_typenoise:movement_typeL5 -0.357 0.212 49.208 -1.685 0.098
Allocationlithium:phasepost3 -0.063 0.100 3373.883 -0.635 0.525
Allocationlithium:phasepost2 0.034 0.099 3375.909 0.338 0.735
Allocationlithium:phasepost1 0.157 0.097 3376.459 1.622 0.105
Allocationlithium:phasepre -0.293 0.096 3379.602 -3.059 0.002
var_typenoise:phasepost3 -0.134 0.099 3373.004 -1.352 0.177
var_typenoise:phasepost2 -0.320 0.099 3374.993 -3.216 0.001
var_typenoise:phasepost1 -0.530 0.099 3382.936 -5.380 0.000
var_typenoise:phasepre -1.242 0.096 3381.821 -12.929 0.000
movement_typeL5:phasepost3 -0.091 0.099 3372.926 -0.923 0.356
movement_typeL5:phasepost2 -0.152 0.099 3374.845 -1.527 0.127
movement_typeL5:phasepost1 -0.063 0.099 3382.630 -0.643 0.520
movement_typeL5:phasepre -0.204 0.096 3381.497 -2.127 0.034
Allocationlithium:var_typenoise:movement_typeL5 -0.025 0.287 51.668 -0.087 0.931
Allocationlithium:var_typenoise:phasepost3 0.110 0.141 3377.044 0.778 0.437
Allocationlithium:var_typenoise:phasepost2 0.210 0.140 3380.014 1.502 0.133
Allocationlithium:var_typenoise:phasepost1 0.060 0.136 3380.058 0.440 0.660
Allocationlithium:var_typenoise:phasepre 0.374 0.135 3385.197 2.777 0.006
Allocationlithium:movement_typeL5:phasepost3 -0.108 0.141 3376.824 -0.765 0.444
Allocationlithium:movement_typeL5:phasepost2 -0.114 0.140 3379.809 -0.817 0.414
Allocationlithium:movement_typeL5:phasepost1 -0.495 0.136 3379.970 -3.626 0.000
Allocationlithium:movement_typeL5:phasepre -0.209 0.135 3384.960 -1.550 0.121
var_typenoise:movement_typeL5:phasepost3 0.115 0.140 3373.786 0.823 0.411
var_typenoise:movement_typeL5:phasepost2 0.361 0.141 3376.481 2.566 0.010
var_typenoise:movement_typeL5:phasepost1 0.284 0.139 3386.439 2.039 0.042
var_typenoise:movement_typeL5:phasepre 0.468 0.136 3385.030 3.447 0.001
Allocationlithium:var_typenoise:movement_typeL5:phasepost3 0.117 0.200 3378.983 0.584 0.559
Allocationlithium:var_typenoise:movement_typeL5:phasepost2 -0.132 0.198 3382.997 -0.665 0.506
Allocationlithium:var_typenoise:movement_typeL5:phasepost1 0.260 0.193 3383.127 1.348 0.178
Allocationlithium:var_typenoise:movement_typeL5:phasepre 0.039 0.191 3389.407 0.203 0.839
df_hlm_lmer_movement_variability$phase <- factor(df_hlm_lmer_movement_variability$phase,
                                          levels = c("pre","post1","post2","post3","post4"))

plot_model(lmer(model_formula,data = df_hlm_lmer_movement_variability),
           type = "pred",terms = c("phase","Allocation"))

df_hlm_lmer_movement_variability$var_type <- relevel(df_hlm_lmer_movement_variability$var_type, ref = "volatility")
df_hlm_lmer_movement_variability$movement_type <- relevel(df_hlm_lmer_movement_variability$movement_type, ref = "M10")
ee2 <- as.data.frame(Effect(c("phase", "Allocation","movement_type","var_type"),mod_movement_var_vmu_M10))

stat.test_s_M10 <- as.data.frame(round(summary(mod_movement_var_s_M10)$coefficients,3))[16:19,]
stat.test_s_L5 <- as.data.frame(round(summary(mod_movement_var_s_L5)$coefficients,3))[16:19,]
stat.test_vmu_M10 <- as.data.frame(round(summary(mod_movement_var_vmu_M10)$coefficients,3))[16:19,]
stat.test_vmu_L5 <- as.data.frame(round(summary(mod_movement_var_vmu_L5)$coefficients,3))[16:19,]

stat.test <- Reduce(rbind,list(stat.test_s_M10,stat.test_s_L5,stat.test_vmu_M10,stat.test_vmu_L5))%>%rename(p_val = 'Pr(>|t|)')
stat.test[c(1:4,9:12),"movement_type"] <- "M10" 
stat.test[c(5:8,13:16),"movement_type"] <- "L5" 

stat.test[1:8,"var_type"] <- "noise" 
stat.test[9:16,"var_type"] <- "volatility" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee2[ee2$phase != "pre" & ee2$Allocation == "lithium",]%>%
  rename(group1 = phase), by = c("movement_type","var_type","group1"))%>%
  mutate(y.position = if_else(var_type == "volatility",ee2$upper[1],if_else(var_type == "noise", ee2$upper[21],0)))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))

#plot noise first
p11a <- ggplot(ee2%>%filter(var_type == "noise"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~movement_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement noise")+
  scale_y_continuous(limits = c(-2,-1),breaks = seq(-2,-1, by = 0.2))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "noise"),label = "p_sig",vjust = -3)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))


#plot volatility first
p11b <- ggplot(ee2%>%filter(var_type == "volatility"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~movement_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement volatility")+
  scale_y_continuous(limits = c(-5.8,-3.8),breaks = seq(-5.8,-3.8, by = 0.4))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

p11a/p11b

ggsave("../figs/hlm_movement_variability.png",width = 7 ,height=7,dpi = 300)

p11b_M10 <- ggplot(ee2%>%filter(var_type == "volatility", movement_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime movement volatility",x="")+
  scale_y_continuous(limits = c(-5.8,-3.8),breaks = seq(-5.8,-3.8, by = 0.4))+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility", movement_type == "M10"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

time

raw value

model_formula_mean <- val ~ Allocation * time_type * phase + posLog + pos_vmu + Age + Gender + (time_type|ID)

df_hlm_lmer_time_mean <- df_hlm%>%
  select(ID,Allocation,phase,day,time_L5_mean,time_M10_mean,Age,Gender,posLog,negLog, pos_vmu)%>%
  pivot_longer(c(time_L5_mean,time_M10_mean),names_pattern =  "time_(.*)_mean",names_to = c("time_type"),values_to = "val")

df_hlm_lmer_time_mean$time_type <- factor(df_hlm_lmer_time_mean$time_type)
df_hlm_lmer_time_mean$time_type <- relevel(df_hlm_lmer_time_mean$time_type, ref = "M10")

mod_time_mean_M10 <- lmer(model_formula_mean,data = df_hlm_lmer_time_mean)
kable(round(summary(mod_time_mean_M10)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 14.566 1.340 100.110 10.869 0.000
Allocationlithium 0.613 0.746 52.525 0.822 0.415
time_typeL5 13.661 0.451 192.100 30.298 0.000
phasepost1 -0.756 0.460 1086.546 -1.645 0.100
phasepost2 -0.468 0.494 1099.553 -0.946 0.344
phasepost3 -0.018 0.455 1098.524 -0.039 0.969
phasepost4 1.022 0.505 1095.268 2.022 0.043
posLog -0.477 0.307 989.204 -1.557 0.120
pos_vmu 0.036 0.093 996.881 0.384 0.701
Age -0.038 0.024 29.776 -1.599 0.120
GenderM 0.058 0.567 31.615 0.103 0.919
Allocationlithium:time_typeL5 -1.164 0.637 218.428 -1.829 0.069
Allocationlithium:phasepost1 -0.504 0.675 1095.374 -0.747 0.455
Allocationlithium:phasepost2 -1.230 0.682 1097.438 -1.803 0.072
Allocationlithium:phasepost3 -1.547 0.658 1098.210 -2.350 0.019
Allocationlithium:phasepost4 -1.691 0.713 1097.265 -2.370 0.018
time_typeL5:phasepost1 0.341 0.642 1091.042 0.532 0.595
time_typeL5:phasepost2 -0.244 0.683 1084.170 -0.358 0.721
time_typeL5:phasepost3 -0.612 0.629 1076.765 -0.974 0.330
time_typeL5:phasepost4 -1.503 0.690 1021.463 -2.177 0.030
Allocationlithium:time_typeL5:phasepost1 0.949 0.943 1093.503 1.007 0.314
Allocationlithium:time_typeL5:phasepost2 1.569 0.948 1083.099 1.655 0.098
Allocationlithium:time_typeL5:phasepost3 1.918 0.907 1059.946 2.116 0.035
Allocationlithium:time_typeL5:phasepost4 1.520 0.973 1041.081 1.562 0.119
df_hlm_lmer_time_mean$time_type <- factor(df_hlm_lmer_time_mean$time_type)
df_hlm_lmer_time_mean$time_type <- relevel(df_hlm_lmer_time_mean$time_type, ref = "L5")

mod_time_mean_L5 <- lmer(model_formula_mean,data = df_hlm_lmer_time_mean)
kable(round(summary(mod_time_mean_L5)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 28.228 1.321 96.001 21.366 0.000
Allocationlithium -0.551 0.678 58.541 -0.812 0.420
time_typeM10 -13.661 0.451 192.105 -30.298 0.000
phasepost1 -0.414 0.459 1088.401 -0.903 0.367
phasepost2 -0.712 0.492 1099.492 -1.446 0.148
phasepost3 -0.630 0.453 1098.461 -1.390 0.165
phasepost4 -0.482 0.501 1079.869 -0.961 0.337
posLog -0.477 0.307 989.206 -1.557 0.120
pos_vmu 0.036 0.093 996.883 0.384 0.701
Age -0.038 0.024 29.777 -1.599 0.120
GenderM 0.058 0.567 31.616 0.103 0.919
Allocationlithium:time_typeM10 1.164 0.637 218.433 1.829 0.069
Allocationlithium:phasepost1 0.445 0.673 1097.623 0.662 0.508
Allocationlithium:phasepost2 0.339 0.680 1098.723 0.499 0.618
Allocationlithium:phasepost3 0.371 0.655 1096.811 0.567 0.571
Allocationlithium:phasepost4 -0.171 0.708 1084.489 -0.241 0.809
time_typeM10:phasepost1 -0.341 0.642 1091.042 -0.532 0.595
time_typeM10:phasepost2 0.244 0.683 1084.170 0.358 0.721
time_typeM10:phasepost3 0.612 0.629 1076.764 0.974 0.330
time_typeM10:phasepost4 1.503 0.690 1021.464 2.177 0.030
Allocationlithium:time_typeM10:phasepost1 -0.949 0.943 1093.503 -1.007 0.314
Allocationlithium:time_typeM10:phasepost2 -1.569 0.948 1083.098 -1.655 0.098
Allocationlithium:time_typeM10:phasepost3 -1.918 0.907 1059.946 -2.116 0.035
Allocationlithium:time_typeM10:phasepost4 -1.520 0.973 1041.081 -1.562 0.119
plot_model(mod_time_mean_M10, type = "pred", terms = c("phase", "Allocation","time_type"))+
  scale_color_manual(values=orangeNavyPalette)+
  theme_bw()+
  labs(title = "predicted value of onset time",y="onset time")
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

df_hlm_lmer_time_mean$time_type <- relevel(df_hlm_lmer_time_mean$time_type, ref = "M10")
ee3 <- as.data.frame(Effect(c("phase", "Allocation","time_type"),mod_time_mean_M10))

stat.test_M10 <- as.data.frame(round(summary(mod_time_mean_M10)$coefficients,3))[13:16,]
stat.test_L5 <- as.data.frame(round(summary(mod_time_mean_L5)$coefficients,3))[13:16,]
stat.test <- rbind(stat.test_L5,stat.test_M10)%>%rename(p_val = 'Pr(>|t|)')
stat.test[1:4,"time_type"] <- "L5" 
stat.test[5:8,"time_type"] <- "M10" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee3[ee3$phase != "pre" & ee3$Allocation == "control",]%>%
  rename(group1 = phase))%>%
  mutate(y.position = if_else(time_type == "L5",ee3$upper[11],if_else(time_type == "M10", ee3$upper[1],0)))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))
## Joining, by = c("time_type", "group1")
p4b <- ggplot(ee3, aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_grid(~time_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="movement onset time (hrs)")+ 
  scale_y_continuous(limits = c(10,30),breaks = seq(10,30, by = 2.5))+
  stat_pvalue_manual(stat.test, label = "p_sig",vjust = -2)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4b

ggsave("../figs/hlm_time_val.png",width = 6 ,height=4,dpi = 300)

# p4a/p4b
# ggsave("../figs/fig2_vals.png",width = 7 ,height=7,dpi = 300)

p4b_M10 <- ggplot(ee3%>%filter(time_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime movement onset time",x="")+ 
  scale_y_continuous(limits = c(8,16),breaks = seq(8, 16, by = 2))+
  stat_pvalue_manual(stat.test%>%filter(time_type == "M10"), label = "p_sig",vjust = -4)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))
p4b_M10

variability

model_formula <- val ~ Allocation * var_type * time_type * phase + posLog + pos_vmu + Age + Gender + (var_type * time_type|ID)

df_hlm_lmer_time_variability <- df_hlm%>%
  select(ID,Allocation,day,phase,time_L5_s,time_L5_vmu, time_M10_s,time_M10_vmu,Age,Gender,posLog, negLog, pos_vmu)%>%
  pivot_longer(c(time_L5_s,time_L5_vmu, time_M10_s,time_M10_vmu),names_pattern =  "time_(.*)_(.*)",names_to = c("time_type","var_type"),values_to = "val")%>%
  mutate(var_type = if_else(var_type == "vmu","volatility","noise"))

df_hlm_lmer_time_variability$var_type <- factor(df_hlm_lmer_time_variability$var_type)
df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "noise")
df_hlm_lmer_time_variability$time_type <- factor(df_hlm_lmer_time_variability$time_type)
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "M10")

mod_time_var_s_M10 <- lmer(model_formula,data = df_hlm_lmer_time_variability)
kable(round(summary(mod_time_var_s_M10)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.274 0.128 95.836 -9.954 0.000
Allocationlithium 0.128 0.093 58.127 1.377 0.174
var_typevolatility -2.649 0.181 37.026 -14.605 0.000
time_typeL5 -0.155 0.103 52.298 -1.497 0.140
phasepost1 0.044 0.058 3382.034 0.747 0.455
phasepost2 0.052 0.060 3397.474 0.864 0.387
phasepost3 0.122 0.060 3393.539 2.047 0.041
phasepost4 0.138 0.061 3392.187 2.268 0.023
posLog 0.068 0.025 2699.760 2.704 0.007
pos_vmu 0.002 0.008 2335.667 0.198 0.843
Age -0.005 0.002 28.986 -2.058 0.049
GenderM 0.019 0.054 31.386 0.358 0.723
Allocationlithium:var_typevolatility -0.203 0.244 37.643 -0.834 0.409
Allocationlithium:time_typeL5 -0.110 0.140 54.807 -0.787 0.435
var_typevolatility:time_typeL5 -0.115 0.174 44.472 -0.660 0.513
Allocationlithium:phasepost1 -0.035 0.080 3384.337 -0.442 0.659
Allocationlithium:phasepost2 -0.092 0.082 3388.245 -1.118 0.264
Allocationlithium:phasepost3 -0.113 0.084 3392.596 -1.343 0.179
Allocationlithium:phasepost4 -0.048 0.085 3396.764 -0.565 0.572
var_typevolatility:phasepost1 -0.661 0.082 3373.348 -8.024 0.000
var_typevolatility:phasepost2 -0.962 0.084 3383.992 -11.407 0.000
var_typevolatility:phasepost3 -1.398 0.084 3381.797 -16.621 0.000
var_typevolatility:phasepost4 -1.557 0.086 3381.627 -18.142 0.000
time_typeL5:phasepost1 -0.251 0.082 3376.343 -3.045 0.002
time_typeL5:phasepost2 -0.237 0.084 3392.179 -2.811 0.005
time_typeL5:phasepost3 -0.367 0.084 3389.237 -4.369 0.000
time_typeL5:phasepost4 -0.408 0.086 3389.084 -4.754 0.000
Allocationlithium:var_typevolatility:time_typeL5 0.186 0.235 46.035 0.792 0.432
Allocationlithium:var_typevolatility:phasepost1 0.171 0.113 3377.229 1.514 0.130
Allocationlithium:var_typevolatility:phasepost2 0.237 0.117 3379.539 2.030 0.042
Allocationlithium:var_typevolatility:phasepost3 0.458 0.119 3381.999 3.858 0.000
Allocationlithium:var_typevolatility:phasepost4 0.381 0.120 3385.310 3.164 0.002
Allocationlithium:time_typeL5:phasepost1 0.036 0.113 3382.715 0.320 0.749
Allocationlithium:time_typeL5:phasepost2 -0.028 0.117 3386.073 -0.244 0.807
Allocationlithium:time_typeL5:phasepost3 0.061 0.118 3389.866 0.512 0.608
Allocationlithium:time_typeL5:phasepost4 0.009 0.120 3394.365 0.073 0.942
var_typevolatility:time_typeL5:phasepost1 0.032 0.116 3374.340 0.274 0.784
var_typevolatility:time_typeL5:phasepost2 0.110 0.119 3387.403 0.926 0.354
var_typevolatility:time_typeL5:phasepost3 0.696 0.119 3384.923 5.857 0.000
var_typevolatility:time_typeL5:phasepost4 0.914 0.121 3384.781 7.529 0.000
Allocationlithium:var_typevolatility:time_typeL5:phasepost1 0.033 0.160 3379.339 0.204 0.838
Allocationlithium:var_typevolatility:time_typeL5:phasepost2 0.257 0.165 3382.111 1.558 0.119
Allocationlithium:var_typevolatility:time_typeL5:phasepost3 -0.377 0.168 3385.311 -2.247 0.025
Allocationlithium:var_typevolatility:time_typeL5:phasepost4 -0.269 0.170 3389.410 -1.579 0.114
df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "volatility")
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "M10")

#same model, different reference level
mod_time_var_vmu_M10 <- lmer(model_formula,data = df_hlm_lmer_time_variability)
kable(round(summary(mod_time_var_vmu_M10)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -3.923 0.187 69.874 -21.007 0.000
Allocationlithium -0.076 0.204 36.120 -0.371 0.713
var_typenoise 2.649 0.181 37.057 14.609 0.000
time_typeL5 -0.270 0.138 41.058 -1.954 0.058
phasepost1 -0.617 0.058 3376.304 -10.553 0.000
phasepost2 -0.910 0.060 3383.686 -15.168 0.000
phasepost3 -1.276 0.060 3380.657 -21.348 0.000
phasepost4 -1.420 0.061 3379.188 -23.302 0.000
posLog 0.068 0.025 2699.634 2.704 0.007
pos_vmu 0.002 0.008 2335.506 0.198 0.843
Age -0.005 0.002 28.990 -2.058 0.049
GenderM 0.019 0.054 31.390 0.357 0.723
Allocationlithium:var_typenoise 0.203 0.244 37.674 0.834 0.409
Allocationlithium:time_typeL5 0.076 0.186 42.229 0.408 0.685
var_typenoise:time_typeL5 0.115 0.174 44.494 0.660 0.513
Allocationlithium:phasepost1 0.136 0.080 3373.991 1.699 0.089
Allocationlithium:phasepost2 0.145 0.083 3375.669 1.751 0.080
Allocationlithium:phasepost3 0.345 0.084 3378.680 4.099 0.000
Allocationlithium:phasepost4 0.333 0.086 3381.195 3.893 0.000
var_typenoise:phasepost1 0.661 0.082 3373.330 8.024 0.000
var_typenoise:phasepost2 0.962 0.084 3383.974 11.407 0.000
var_typenoise:phasepost3 1.398 0.084 3381.779 16.621 0.000
var_typenoise:phasepost4 1.557 0.086 3381.610 18.141 0.000
time_typeL5:phasepost1 -0.219 0.082 3372.806 -2.657 0.008
time_typeL5:phasepost2 -0.126 0.084 3383.055 -1.498 0.134
time_typeL5:phasepost3 0.329 0.084 3381.127 3.914 0.000
time_typeL5:phasepost4 0.506 0.086 3380.988 5.893 0.000
Allocationlithium:var_typenoise:time_typeL5 -0.186 0.235 46.058 -0.792 0.432
Allocationlithium:var_typenoise:phasepost1 -0.171 0.113 3377.211 -1.514 0.130
Allocationlithium:var_typenoise:phasepost2 -0.237 0.117 3379.521 -2.030 0.042
Allocationlithium:var_typenoise:phasepost3 -0.458 0.119 3381.982 -3.858 0.000
Allocationlithium:var_typenoise:phasepost4 -0.381 0.120 3385.293 -3.164 0.002
Allocationlithium:time_typeL5:phasepost1 0.069 0.113 3376.569 0.608 0.543
Allocationlithium:time_typeL5:phasepost2 0.228 0.117 3378.776 1.959 0.050
Allocationlithium:time_typeL5:phasepost3 -0.316 0.119 3381.270 -2.664 0.008
Allocationlithium:time_typeL5:phasepost4 -0.260 0.120 3384.668 -2.159 0.031
var_typenoise:time_typeL5:phasepost1 -0.032 0.116 3374.322 -0.274 0.784
var_typenoise:time_typeL5:phasepost2 -0.110 0.119 3387.385 -0.926 0.354
var_typenoise:time_typeL5:phasepost3 -0.696 0.119 3384.905 -5.857 0.000
var_typenoise:time_typeL5:phasepost4 -0.914 0.121 3384.763 -7.529 0.000
Allocationlithium:var_typenoise:time_typeL5:phasepost1 -0.033 0.160 3379.321 -0.204 0.838
Allocationlithium:var_typenoise:time_typeL5:phasepost2 -0.257 0.165 3382.093 -1.558 0.119
Allocationlithium:var_typenoise:time_typeL5:phasepost3 0.377 0.168 3385.293 2.247 0.025
Allocationlithium:var_typenoise:time_typeL5:phasepost4 0.269 0.170 3389.392 1.579 0.114
df_hlm_lmer_time_variability$Allocation <- factor(df_hlm_lmer_time_variability$Allocation, levels = c("lithium","control"))
kable(round(summary(lmer(model_formula,data = df_hlm_lmer_time_variability))$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -3.999 0.169 76.963 -23.623 0.000
Allocationcontrol 0.076 0.204 36.098 0.370 0.713
var_typenoise 2.852 0.163 38.441 17.534 0.000
time_typeL5 -0.194 0.125 43.712 -1.554 0.127
phasepost1 -0.481 0.055 3379.372 -8.729 0.000
phasepost2 -0.766 0.057 3376.181 -13.391 0.000
phasepost3 -0.931 0.059 3377.152 -15.724 0.000
phasepost4 -1.087 0.060 3382.310 -18.151 0.000
posLog 0.068 0.025 2699.721 2.704 0.007
pos_vmu 0.002 0.008 2335.612 0.198 0.843
Age -0.005 0.002 28.989 -2.058 0.049
GenderM 0.019 0.054 31.389 0.358 0.723
Allocationcontrol:var_typenoise -0.203 0.244 37.652 -0.834 0.409
Allocationcontrol:time_typeL5 -0.076 0.186 42.207 -0.408 0.685
var_typenoise:time_typeL5 -0.071 0.158 48.058 -0.452 0.653
Allocationcontrol:phasepost1 -0.136 0.080 3374.004 -1.699 0.089
Allocationcontrol:phasepost2 -0.145 0.083 3375.682 -1.751 0.080
Allocationcontrol:phasepost3 -0.345 0.084 3378.693 -4.099 0.000
Allocationcontrol:phasepost4 -0.333 0.086 3381.208 -3.893 0.000
var_typenoise:phasepost1 0.490 0.078 3381.263 6.303 0.000
var_typenoise:phasepost2 0.725 0.081 3374.115 9.010 0.000
var_typenoise:phasepost3 0.940 0.084 3382.479 11.241 0.000
var_typenoise:phasepost4 1.176 0.084 3389.126 13.924 0.000
time_typeL5:phasepost1 -0.150 0.078 3380.515 -1.932 0.053
time_typeL5:phasepost2 0.102 0.081 3373.612 1.267 0.205
time_typeL5:phasepost3 0.013 0.084 3381.464 0.158 0.874
time_typeL5:phasepost4 0.246 0.085 3388.174 2.910 0.004
Allocationcontrol:var_typenoise:time_typeL5 0.186 0.235 46.042 0.792 0.432
Allocationcontrol:var_typenoise:phasepost1 0.171 0.113 3377.226 1.514 0.130
Allocationcontrol:var_typenoise:phasepost2 0.237 0.117 3379.535 2.030 0.042
Allocationcontrol:var_typenoise:phasepost3 0.458 0.119 3381.996 3.858 0.000
Allocationcontrol:var_typenoise:phasepost4 0.381 0.120 3385.307 3.164 0.002
Allocationcontrol:time_typeL5:phasepost1 -0.069 0.113 3376.582 -0.608 0.543
Allocationcontrol:time_typeL5:phasepost2 -0.228 0.117 3378.789 -1.959 0.050
Allocationcontrol:time_typeL5:phasepost3 0.316 0.119 3381.282 2.664 0.008
Allocationcontrol:time_typeL5:phasepost4 0.260 0.120 3384.680 2.159 0.031
var_typenoise:time_typeL5:phasepost1 -0.065 0.110 3384.360 -0.588 0.557
var_typenoise:time_typeL5:phasepost2 -0.367 0.114 3375.348 -3.227 0.001
var_typenoise:time_typeL5:phasepost3 -0.320 0.118 3385.693 -2.703 0.007
var_typenoise:time_typeL5:phasepost4 -0.645 0.119 3393.402 -5.401 0.000
Allocationcontrol:var_typenoise:time_typeL5:phasepost1 0.033 0.160 3379.334 0.204 0.838
Allocationcontrol:var_typenoise:time_typeL5:phasepost2 0.257 0.165 3382.107 1.558 0.119
Allocationcontrol:var_typenoise:time_typeL5:phasepost3 -0.377 0.168 3385.306 -2.247 0.025
Allocationcontrol:var_typenoise:time_typeL5:phasepost4 -0.269 0.170 3389.406 -1.579 0.114
df_hlm_lmer_time_variability$Allocation <- relevel(df_hlm_lmer_time_variability$Allocation, ref = "control")

# #same model, different reference level
df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "noise")
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "L5")

mod_time_var_s_L5 <- lmer(model_formula,data = df_hlm_lmer_time_variability)
kable(round(summary(mod_time_var_s_L5)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.429 0.133 99.677 -10.773 0.000
Allocationlithium 0.017 0.104 50.012 0.167 0.868
var_typevolatility -2.764 0.146 40.036 -18.941 0.000
time_typeM10 0.155 0.103 52.278 1.497 0.140
phasepost1 -0.207 0.058 3380.288 -3.542 0.000
phasepost2 -0.185 0.060 3395.223 -3.089 0.002
phasepost3 -0.245 0.060 3390.979 -4.104 0.000
phasepost4 -0.270 0.061 3389.483 -4.435 0.000
posLog 0.068 0.025 2699.736 2.703 0.007
pos_vmu 0.002 0.008 2335.626 0.198 0.843
Age -0.005 0.002 28.991 -2.058 0.049
GenderM 0.019 0.054 31.391 0.358 0.723
Allocationlithium:var_typevolatility -0.017 0.196 41.048 -0.087 0.931
Allocationlithium:time_typeM10 0.110 0.140 54.786 0.787 0.435
var_typevolatility:time_typeM10 0.115 0.174 44.465 0.659 0.513
Allocationlithium:phasepost1 0.001 0.080 3381.506 0.011 0.992
Allocationlithium:phasepost2 -0.121 0.082 3384.946 -1.463 0.144
Allocationlithium:phasepost3 -0.052 0.084 3389.551 -0.619 0.536
Allocationlithium:phasepost4 -0.039 0.085 3393.976 -0.461 0.645
var_typevolatility:phasepost1 -0.629 0.082 3374.633 -7.637 0.000
var_typevolatility:phasepost2 -0.852 0.084 3388.044 -10.101 0.000
var_typevolatility:phasepost3 -0.701 0.084 3385.452 -8.344 0.000
var_typevolatility:phasepost4 -0.644 0.086 3385.276 -7.500 0.000
time_typeM10:phasepost1 0.251 0.082 3376.346 3.045 0.002
time_typeM10:phasepost2 0.237 0.084 3392.181 2.811 0.005
time_typeM10:phasepost3 0.367 0.084 3389.240 4.369 0.000
time_typeM10:phasepost4 0.408 0.086 3389.087 4.755 0.000
Allocationlithium:var_typevolatility:time_typeM10 -0.186 0.235 46.028 -0.792 0.432
Allocationlithium:var_typevolatility:phasepost1 0.204 0.113 3379.748 1.803 0.071
Allocationlithium:var_typevolatility:phasepost2 0.494 0.117 3382.631 4.234 0.000
Allocationlithium:var_typevolatility:phasepost3 0.081 0.119 3385.807 0.682 0.495
Allocationlithium:var_typevolatility:phasepost4 0.112 0.120 3389.898 0.931 0.352
Allocationlithium:time_typeM10:phasepost1 -0.036 0.113 3382.718 -0.320 0.749
Allocationlithium:time_typeM10:phasepost2 0.028 0.117 3386.075 0.244 0.807
Allocationlithium:time_typeM10:phasepost3 -0.061 0.118 3389.869 -0.512 0.608
Allocationlithium:time_typeM10:phasepost4 -0.009 0.120 3394.367 -0.073 0.942
var_typevolatility:time_typeM10:phasepost1 -0.032 0.116 3374.343 -0.274 0.784
var_typevolatility:time_typeM10:phasepost2 -0.110 0.119 3387.406 -0.926 0.354
var_typevolatility:time_typeM10:phasepost3 -0.696 0.119 3384.926 -5.857 0.000
var_typevolatility:time_typeM10:phasepost4 -0.914 0.121 3384.784 -7.529 0.000
Allocationlithium:var_typevolatility:time_typeM10:phasepost1 -0.033 0.160 3379.342 -0.204 0.838
Allocationlithium:var_typevolatility:time_typeM10:phasepost2 -0.257 0.165 3382.114 -1.558 0.119
Allocationlithium:var_typevolatility:time_typeM10:phasepost3 0.377 0.168 3385.314 2.247 0.025
Allocationlithium:var_typevolatility:time_typeM10:phasepost4 0.269 0.170 3389.413 1.579 0.114
df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "volatility")

df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "L5")
mod_time_var_vmu_L5 <- lmer(model_formula,data = df_hlm_lmer_time_variability)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00215358 (tol = 0.002, component 1)
kable(round(summary(mod_time_var_vmu_L5)$coefficients,3))
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -4.193 0.183 71.618 -22.874 0.000
Allocationlithium 0.000 0.199 36.374 0.001 0.999
var_typenoise 2.764 0.146 40.064 18.945 0.000
time_typeM10 0.270 0.138 41.047 1.953 0.058
phasepost1 -0.836 0.058 3376.356 -14.295 0.000
phasepost2 -1.037 0.060 3384.210 -17.275 0.000
phasepost3 -0.946 0.060 3381.137 -15.839 0.000
phasepost4 -0.914 0.061 3379.763 -14.997 0.000
posLog 0.068 0.025 2699.711 2.703 0.007
pos_vmu 0.002 0.008 2335.602 0.198 0.843
Age -0.005 0.002 28.991 -2.058 0.049
GenderM 0.019 0.054 31.391 0.358 0.723
Allocationlithium:var_typenoise 0.017 0.196 41.077 0.087 0.931
Allocationlithium:time_typeM10 -0.076 0.186 42.218 -0.408 0.685
var_typenoise:time_typeM10 -0.115 0.174 44.490 -0.660 0.513
Allocationlithium:phasepost1 0.205 0.080 3374.255 2.559 0.011
Allocationlithium:phasepost2 0.373 0.083 3375.930 4.519 0.000
Allocationlithium:phasepost3 0.029 0.084 3379.156 0.343 0.732
Allocationlithium:phasepost4 0.073 0.085 3381.853 0.851 0.395
var_typenoise:phasepost1 0.629 0.082 3374.619 7.637 0.000
var_typenoise:phasepost2 0.852 0.084 3388.031 10.101 0.000
var_typenoise:phasepost3 0.701 0.084 3385.439 8.344 0.000
var_typenoise:phasepost4 0.644 0.086 3385.263 7.500 0.000
time_typeM10:phasepost1 0.219 0.082 3372.812 2.657 0.008
time_typeM10:phasepost2 0.126 0.084 3383.060 1.498 0.134
time_typeM10:phasepost3 -0.329 0.084 3381.133 -3.914 0.000
time_typeM10:phasepost4 -0.506 0.086 3380.994 -5.893 0.000
Allocationlithium:var_typenoise:time_typeM10 0.186 0.235 46.054 0.792 0.432
Allocationlithium:var_typenoise:phasepost1 -0.204 0.113 3379.735 -1.803 0.071
Allocationlithium:var_typenoise:phasepost2 -0.494 0.117 3382.618 -4.234 0.000
Allocationlithium:var_typenoise:phasepost3 -0.081 0.119 3385.794 -0.682 0.495
Allocationlithium:var_typenoise:phasepost4 -0.112 0.120 3389.886 -0.931 0.352
Allocationlithium:time_typeM10:phasepost1 -0.069 0.113 3376.574 -0.608 0.543
Allocationlithium:time_typeM10:phasepost2 -0.228 0.117 3378.782 -1.959 0.050
Allocationlithium:time_typeM10:phasepost3 0.316 0.119 3381.275 2.664 0.008
Allocationlithium:time_typeM10:phasepost4 0.260 0.120 3384.673 2.159 0.031
var_typenoise:time_typeM10:phasepost1 0.032 0.116 3374.329 0.274 0.784
var_typenoise:time_typeM10:phasepost2 0.110 0.119 3387.394 0.926 0.354
var_typenoise:time_typeM10:phasepost3 0.696 0.119 3384.914 5.857 0.000
var_typenoise:time_typeM10:phasepost4 0.914 0.121 3384.771 7.529 0.000
Allocationlithium:var_typenoise:time_typeM10:phasepost1 0.033 0.160 3379.328 0.204 0.838
Allocationlithium:var_typenoise:time_typeM10:phasepost2 0.257 0.165 3382.101 1.558 0.119
Allocationlithium:var_typenoise:time_typeM10:phasepost3 -0.377 0.168 3385.301 -2.247 0.025
Allocationlithium:var_typenoise:time_typeM10:phasepost4 -0.269 0.170 3389.402 -1.579 0.114
plot_model(mod_time_var_s_M10, type = "pred", terms = c("phase", "Allocation","time_type","var_type"))+
  labs(title="predicted value of onset time variability",y="noise")

ggsave("../figs/hlm_time_variability_comparison.png",width = 7 ,height=7,dpi = 300)
df_hlm_lmer_time_variability$var_type <- relevel(df_hlm_lmer_time_variability$var_type, ref = "volatility")
df_hlm_lmer_time_variability$time_type <- relevel(df_hlm_lmer_time_variability$time_type, ref = "M10")
ee4 <- as.data.frame(Effect(c("phase", "Allocation","time_type","var_type"),mod_time_var_vmu_M10))

stat.test_s_M10 <- as.data.frame(round(summary(mod_time_var_s_M10)$coefficients,3))[16:19,]
stat.test_s_L5 <- as.data.frame(round(summary(mod_time_var_s_L5)$coefficients,3))[16:19,]
stat.test_vmu_M10 <- as.data.frame(round(summary(mod_time_var_vmu_M10)$coefficients,3))[16:19,]
stat.test_vmu_L5 <- as.data.frame(round(summary(mod_time_var_vmu_L5)$coefficients,3))[16:19,]

stat.test <- Reduce(rbind,list(stat.test_s_M10,stat.test_s_L5,stat.test_vmu_M10,stat.test_vmu_L5))%>%rename(p_val = 'Pr(>|t|)')
stat.test[c(1:4,9:12),"time_type"] <- "M10" 
stat.test[c(5:8,13:16),"time_type"] <- "L5" 

stat.test[1:8,"var_type"] <- "noise" 
stat.test[9:16,"var_type"] <- "volatility" 

stat.test$group1 <- c("post1","post2","post3","post4")
stat.test$group2 <- c("post1","post2","post3","post4")

stat.test <- left_join(stat.test,ee4[ee4$phase != "pre" & ee4$Allocation == "lithium",]%>%
  rename(group1 = phase), by = c("time_type","var_type","group1"))%>%
  mutate(y.position = if_else(var_type == "volatility" & time_type == "L5",ee4$upper[11],if_else(var_type == "volatility" & time_type == "M10", ee4$upper[1],if_else(var_type == "noise" & time_type == "L5",ee4$upper[36],if_else(var_type == "noise" & time_type == "M10",ee4$upper[21],0)))))%>%
  mutate(p_sig = if_else(p_val < 0.001,"***", if_else(p_val < 0.01, "**", if_else(p_val < 0.05, "*", "ns")) ))

#plot noise first
p12a <- ggplot(ee4%>%filter(var_type == "noise"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~time_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement noise")+
  stat_pvalue_manual(stat.test%>%filter(var_type == "noise"),label = "p_sig",vjust = -1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

#plot volatility first
p12b <- ggplot(ee4%>%filter(var_type == "volatility"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  facet_wrap(~time_type)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank())+
  labs(title = "",y="movement volatility")+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

p12a/p12b

ggsave("../figs/hlm_time_variability.png",width = 7 ,height=7,dpi = 300)

p12b_M10 <- ggplot(ee4%>%filter(var_type == "volatility", time_type == "M10"), aes(x = phase, y = fit, color = Allocation))+
  scale_color_manual(values=orangeNavyPalette)+
  geom_errorbar(aes(ymin = lower,ymax=upper),width=0.1,color = "black")+
    geom_point(size = 3)+
  theme(panel.grid.major.y = element_line(color = "black",size = 0.2),
        panel.grid.minor.y = element_blank(),
         panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle=45))+
  labs(title = "",y="daytime onset time volatility",x="")+
  stat_pvalue_manual(stat.test%>%filter(var_type == "volatility", time_type == "M10"),label = "p_sig",vjust = 1)+
  scale_x_discrete(labels = c("run-in","week 1","week 2","week 3","week 4"))

#manuscript graphs

p4a_M10 + p11b_M10+
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/movement.png",width = 7 ,height=3,dpi = 300)

p4b_M10 + p12b_M10+
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/time.png",width = 7 ,height=3,dpi = 300)

(p4a_M10 + p11b_M10) / (p4b_M10 + p12b_M10) +
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/FIG4_movement_n_time.png",width = 7 ,height=6,dpi = 300)
  
(p8 + p9) / (p4a_M10 + p11b_M10) / (p4b_M10 + p12b_M10) +
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/movement_n_time_n_pa.png",width = 7 ,height=10,dpi = 300)

correlating change in vmu and mean

df_hlm$phase <- factor(df_hlm$phase, labels = c("run-in","week 1","week 2","week 3","week 4"))

cA1 <- lmer(time_M10_mean ~ Allocation * movement_M10_mean * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cA1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: time_M10_mean ~ Allocation * movement_M10_mean * phase + Age +  
##     Gender + (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 4133.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3191 -0.7139 -0.0724  0.7163  2.8656 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  2.349   1.533   
##  Residual             11.257   3.355   
## Number of obs: 767, groups:  ID, 34
## 
## Fixed effects:
##                                                   Estimate Std. Error
## (Intercept)                                      14.083237   1.491655
## Allocationlithium                                -0.306840   1.432668
## movement_M10_mean                                -0.014464   0.017957
## phaseweek 1                                      -1.711169   1.547788
## phaseweek 2                                      -2.010896   1.589343
## phaseweek 3                                      -2.235538   1.399565
## phaseweek 4                                      -1.483085   1.616221
## Age                                              -0.040895   0.025603
## GenderM                                           0.803187   0.603491
## Allocationlithium:movement_M10_mean               0.006565   0.021387
## Allocationlithium:phaseweek 1                     1.473474   2.040221
## Allocationlithium:phaseweek 2                     1.561965   1.925902
## Allocationlithium:phaseweek 3                     1.316191   1.811787
## Allocationlithium:phaseweek 4                     0.924409   2.045464
## movement_M10_mean:phaseweek 1                     0.012929   0.025778
## movement_M10_mean:phaseweek 2                     0.026941   0.024114
## movement_M10_mean:phaseweek 3                     0.031050   0.021608
## movement_M10_mean:phaseweek 4                     0.038700   0.025195
## Allocationlithium:movement_M10_mean:phaseweek 1  -0.028480   0.033281
## Allocationlithium:movement_M10_mean:phaseweek 2  -0.040369   0.030056
## Allocationlithium:movement_M10_mean:phaseweek 3  -0.035999   0.029345
## Allocationlithium:movement_M10_mean:phaseweek 4  -0.035553   0.033116
##                                                         df t value Pr(>|t|)    
## (Intercept)                                     110.389365   9.441 7.25e-16 ***
## Allocationlithium                               387.029151  -0.214    0.831    
## movement_M10_mean                               723.112305  -0.805    0.421    
## phaseweek 1                                     724.774360  -1.106    0.269    
## phaseweek 2                                     721.807503  -1.265    0.206    
## phaseweek 3                                     732.081078  -1.597    0.111    
## phaseweek 4                                     734.518937  -0.918    0.359    
## Age                                              30.765009  -1.597    0.120    
## GenderM                                          29.895186   1.331    0.193    
## Allocationlithium:movement_M10_mean             738.006859   0.307    0.759    
## Allocationlithium:phaseweek 1                   729.009206   0.722    0.470    
## Allocationlithium:phaseweek 2                   722.846139   0.811    0.418    
## Allocationlithium:phaseweek 3                   730.096092   0.726    0.468    
## Allocationlithium:phaseweek 4                   736.339936   0.452    0.651    
## movement_M10_mean:phaseweek 1                   728.339714   0.502    0.616    
## movement_M10_mean:phaseweek 2                   720.840497   1.117    0.264    
## movement_M10_mean:phaseweek 3                   732.617335   1.437    0.151    
## movement_M10_mean:phaseweek 4                   740.436632   1.536    0.125    
## Allocationlithium:movement_M10_mean:phaseweek 1 729.139129  -0.856    0.392    
## Allocationlithium:movement_M10_mean:phaseweek 2 721.536396  -1.343    0.180    
## Allocationlithium:movement_M10_mean:phaseweek 3 729.844951  -1.227    0.220    
## Allocationlithium:movement_M10_mean:phaseweek 4 737.476413  -1.074    0.283    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
cA2 <- lmer(pos_mean ~ Allocation * movement_M10_mean * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cA2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_mean ~ Allocation * movement_M10_mean * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 2997
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8470 -0.5705 -0.0628  0.5891  3.2847 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 8.772    2.962   
##  Residual             8.530    2.921   
## Number of obs: 575, groups:  ID, 34
## 
## Fixed effects:
##                                                   Estimate Std. Error
## (Intercept)                                       9.266784   2.098319
## Allocationlithium                                 3.476211   1.653969
## movement_M10_mean                                 0.038176   0.016858
## phaseweek 1                                      -3.554955   1.499697
## phaseweek 2                                      -0.985395   1.493013
## phaseweek 3                                       1.667830   1.305060
## phaseweek 4                                       2.500531   1.477875
## Age                                               0.024768   0.045479
## GenderM                                          -3.611984   1.078351
## Allocationlithium:movement_M10_mean              -0.034361   0.020190
## Allocationlithium:phaseweek 1                    -0.454923   2.004505
## Allocationlithium:phaseweek 2                    -1.433871   1.882547
## Allocationlithium:phaseweek 3                    -3.343205   1.796356
## Allocationlithium:phaseweek 4                    -3.787107   1.995389
## movement_M10_mean:phaseweek 1                     0.034753   0.025166
## movement_M10_mean:phaseweek 2                     0.005019   0.022474
## movement_M10_mean:phaseweek 3                    -0.023569   0.020040
## movement_M10_mean:phaseweek 4                    -0.043633   0.023100
## Allocationlithium:movement_M10_mean:phaseweek 1   0.002307   0.032607
## Allocationlithium:movement_M10_mean:phaseweek 2   0.005677   0.028800
## Allocationlithium:movement_M10_mean:phaseweek 3   0.031761   0.029160
## Allocationlithium:movement_M10_mean:phaseweek 4   0.039144   0.032336
##                                                         df t value Pr(>|t|)    
## (Intercept)                                      50.282666   4.416 5.34e-05 ***
## Allocationlithium                               139.115978   2.102  0.03738 *  
## movement_M10_mean                               544.946598   2.265  0.02393 *  
## phaseweek 1                                     527.556294  -2.370  0.01813 *  
## phaseweek 2                                     526.220935  -0.660  0.50954    
## phaseweek 3                                     529.156148   1.278  0.20182    
## phaseweek 4                                     530.576186   1.692  0.09124 .  
## Age                                              30.619426   0.545  0.58998    
## GenderM                                          30.873126  -3.350  0.00215 ** 
## Allocationlithium:movement_M10_mean             541.477980  -1.702  0.08935 .  
## Allocationlithium:phaseweek 1                   528.689951  -0.227  0.82055    
## Allocationlithium:phaseweek 2                   526.941055  -0.762  0.44660    
## Allocationlithium:phaseweek 3                   528.071771  -1.861  0.06329 .  
## Allocationlithium:phaseweek 4                   534.255909  -1.898  0.05824 .  
## movement_M10_mean:phaseweek 1                   529.065792   1.381  0.16787    
## movement_M10_mean:phaseweek 2                   525.556403   0.223  0.82338    
## movement_M10_mean:phaseweek 3                   528.844678  -1.176  0.24007    
## movement_M10_mean:phaseweek 4                   534.631607  -1.889  0.05945 .  
## Allocationlithium:movement_M10_mean:phaseweek 1 529.057599   0.071  0.94363    
## Allocationlithium:movement_M10_mean:phaseweek 2 525.966282   0.197  0.84382    
## Allocationlithium:movement_M10_mean:phaseweek 3 527.430213   1.089  0.27657    
## Allocationlithium:movement_M10_mean:phaseweek 4 533.493385   1.211  0.22661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
cA3 <- lmer(pos_mean ~ Allocation * time_M10_mean * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cA3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_mean ~ Allocation * time_M10_mean * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 2969.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.63961 -0.57817 -0.07669  0.61344  2.94435 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 8.965    2.994   
##  Residual             8.688    2.948   
## Number of obs: 575, groups:  ID, 34
## 
## Fixed effects:
##                                              Estimate Std. Error        df
## (Intercept)                                  14.34854    2.27998  66.20570
## Allocationlithium                             0.28204    2.07303 263.32115
## time_M10_mean                                -0.22419    0.10113 538.41941
## phaseweek 1                                  -6.73621    1.82388 527.40012
## phaseweek 2                                  -2.87334    1.91234 527.77995
## phaseweek 3                                  -2.31736    1.80492 532.64005
## phaseweek 4                                  -1.98327    1.99983 529.49126
## Age                                           0.01628    0.04601  30.63413
## GenderM                                      -3.54187    1.08844  30.62932
## Allocationlithium:time_M10_mean               0.11098    0.13822 536.02124
## Allocationlithium:phaseweek 1                 2.75382    2.49878 529.07945
## Allocationlithium:phaseweek 2                -0.37921    2.64428 528.12459
## Allocationlithium:phaseweek 3                 2.09386    2.51152 530.81644
## Allocationlithium:phaseweek 4                -2.09369    2.69668 528.82996
## time_M10_mean:phaseweek 1                     0.45386    0.15449 527.13715
## time_M10_mean:phaseweek 2                     0.21242    0.16199 529.56932
## time_M10_mean:phaseweek 3                     0.23268    0.14976 535.76168
## time_M10_mean:phaseweek 4                     0.18096    0.15957 531.87043
## Allocationlithium:time_M10_mean:phaseweek 1  -0.28448    0.20518 528.96994
## Allocationlithium:time_M10_mean:phaseweek 2  -0.09835    0.22198 530.04149
## Allocationlithium:time_M10_mean:phaseweek 3  -0.33724    0.20386 532.29166
## Allocationlithium:time_M10_mean:phaseweek 4   0.02463    0.21174 529.64031
##                                             t value Pr(>|t|)    
## (Intercept)                                   6.293 2.83e-08 ***
## Allocationlithium                             0.136 0.891882    
## time_M10_mean                                -2.217 0.027055 *  
## phaseweek 1                                  -3.693 0.000244 ***
## phaseweek 2                                  -1.503 0.133559    
## phaseweek 3                                  -1.284 0.199730    
## phaseweek 4                                  -0.992 0.321788    
## Age                                           0.354 0.725791    
## GenderM                                      -3.254 0.002773 ** 
## Allocationlithium:time_M10_mean               0.803 0.422361    
## Allocationlithium:phaseweek 1                 1.102 0.270935    
## Allocationlithium:phaseweek 2                -0.143 0.886024    
## Allocationlithium:phaseweek 3                 0.834 0.404825    
## Allocationlithium:phaseweek 4                -0.776 0.437863    
## time_M10_mean:phaseweek 1                     2.938 0.003451 ** 
## time_M10_mean:phaseweek 2                     1.311 0.190330    
## time_M10_mean:phaseweek 3                     1.554 0.120867    
## time_M10_mean:phaseweek 4                     1.134 0.257285    
## Allocationlithium:time_M10_mean:phaseweek 1  -1.387 0.166172    
## Allocationlithium:time_M10_mean:phaseweek 2  -0.443 0.657925    
## Allocationlithium:time_M10_mean:phaseweek 3  -1.654 0.098663 .  
## Allocationlithium:time_M10_mean:phaseweek 4   0.116 0.907456    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
cB1 <- lmer(time_M10_vmu ~ Allocation * movement_M10_vmu * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cB1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: time_M10_vmu ~ Allocation * movement_M10_vmu * phase + Age +  
##     Gender + (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 1700.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1077 -0.6310 -0.0649  0.5670  3.9057 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.3303   0.5747  
##  Residual             0.1951   0.4417  
## Number of obs: 1243, groups:  ID, 34
## 
## Fixed effects:
##                                                  Estimate Std. Error         df
## (Intercept)                                      -2.17590    0.43103   69.79928
## Allocationlithium                                 0.44302    0.37654  286.96038
## movement_M10_vmu                                  0.32468    0.05884 1196.34241
## phaseweek 1                                       0.03898    0.42094 1196.79933
## phaseweek 2                                       0.03499    0.42611 1197.60737
## phaseweek 3                                      -0.13169    0.46659 1198.93349
## phaseweek 4                                      -0.60940    0.48631 1195.77350
## Age                                              -0.01257    0.00858   29.97599
## GenderM                                           0.07032    0.20314   30.10371
## Allocationlithium:movement_M10_vmu                0.14081    0.07449 1199.04329
## Allocationlithium:phaseweek 1                    -1.25961    0.49930 1195.80379
## Allocationlithium:phaseweek 2                    -1.01360    0.51375 1196.67988
## Allocationlithium:phaseweek 3                    -1.76943    0.54808 1198.21475
## Allocationlithium:phaseweek 4                    -1.29420    0.56350 1196.07011
## movement_M10_vmu:phaseweek 1                      0.06931    0.08793 1195.71729
## movement_M10_vmu:phaseweek 2                      0.11231    0.08711 1196.12480
## movement_M10_vmu:phaseweek 3                      0.13962    0.09317 1197.22463
## movement_M10_vmu:phaseweek 4                      0.06966    0.09559 1194.70554
## Allocationlithium:movement_M10_vmu:phaseweek 1   -0.28495    0.10686 1194.77906
## Allocationlithium:movement_M10_vmu:phaseweek 2   -0.22793    0.10675 1195.21961
## Allocationlithium:movement_M10_vmu:phaseweek 3   -0.43293    0.11161 1196.60859
## Allocationlithium:movement_M10_vmu:phaseweek 4   -0.33428    0.11352 1195.18760
##                                                t value Pr(>|t|)    
## (Intercept)                                     -5.048 3.41e-06 ***
## Allocationlithium                                1.177 0.240346    
## movement_M10_vmu                                 5.518 4.20e-08 ***
## phaseweek 1                                      0.093 0.926234    
## phaseweek 2                                      0.082 0.934577    
## phaseweek 3                                     -0.282 0.777814    
## phaseweek 4                                     -1.253 0.210414    
## Age                                             -1.465 0.153452    
## GenderM                                          0.346 0.731626    
## Allocationlithium:movement_M10_vmu               1.890 0.058965 .  
## Allocationlithium:phaseweek 1                   -2.523 0.011773 *  
## Allocationlithium:phaseweek 2                   -1.973 0.048732 *  
## Allocationlithium:phaseweek 3                   -3.228 0.001279 ** 
## Allocationlithium:phaseweek 4                   -2.297 0.021807 *  
## movement_M10_vmu:phaseweek 1                     0.788 0.430709    
## movement_M10_vmu:phaseweek 2                     1.289 0.197536    
## movement_M10_vmu:phaseweek 3                     1.499 0.134256    
## movement_M10_vmu:phaseweek 4                     0.729 0.466259    
## Allocationlithium:movement_M10_vmu:phaseweek 1  -2.667 0.007765 ** 
## Allocationlithium:movement_M10_vmu:phaseweek 2  -2.135 0.032945 *  
## Allocationlithium:movement_M10_vmu:phaseweek 3  -3.879 0.000111 ***
## Allocationlithium:movement_M10_vmu:phaseweek 4  -2.945 0.003296 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
pC1 <- plot_model(cB1, type = "pred", terms = c("movement_M10_vmu","Allocation","phase"), show.data = F)+
  scale_color_manual(values=orangeNavyPalette)+
  scale_fill_manual(values=orangeNavyPalette)+
  labs(title = "",
       x = "daytime movement volatility",
       y = "movement onset time volatility")
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
pC1$facet$params$nrow=1

cB2 <- lmer(pos_vmu ~ Allocation * movement_M10_vmu * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cB2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_vmu ~ Allocation * movement_M10_vmu * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 3021.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5610 -0.6593 -0.0498  0.6285  4.0737 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.8128   0.9016  
##  Residual             0.5784   0.7605  
## Number of obs: 1243, groups:  ID, 34
## 
## Fixed effects:
##                                                  Estimate Std. Error         df
## (Intercept)                                    -1.462e+00  7.004e-01  7.951e+01
## Allocationlithium                              -1.389e+00  6.316e-01  3.524e+02
## movement_M10_vmu                                5.542e-01  1.013e-01  1.197e+03
## phaseweek 1                                    -2.307e+00  7.245e-01  1.198e+03
## phaseweek 2                                    -4.173e+00  7.334e-01  1.199e+03
## phaseweek 3                                    -4.082e+00  8.030e-01  1.200e+03
## phaseweek 4                                    -1.310e+00  8.371e-01  1.197e+03
## Age                                             5.872e-03  1.348e-02  2.998e+01
## GenderM                                        -6.601e-01  3.193e-01  3.013e+01
## Allocationlithium:movement_M10_vmu             -3.097e-01  1.282e-01  1.200e+03
## Allocationlithium:phaseweek 1                  -5.764e-02  8.594e-01  1.197e+03
## Allocationlithium:phaseweek 2                   7.993e-01  8.843e-01  1.198e+03
## Allocationlithium:phaseweek 3                   3.074e+00  9.433e-01  1.200e+03
## Allocationlithium:phaseweek 4                   1.813e-01  9.699e-01  1.197e+03
## movement_M10_vmu:phaseweek 1                   -4.397e-01  1.514e-01  1.197e+03
## movement_M10_vmu:phaseweek 2                   -7.595e-01  1.499e-01  1.197e+03
## movement_M10_vmu:phaseweek 3                   -7.462e-01  1.604e-01  1.198e+03
## movement_M10_vmu:phaseweek 4                   -2.459e-01  1.645e-01  1.195e+03
## Allocationlithium:movement_M10_vmu:phaseweek 1 -5.375e-03  1.839e-01  1.195e+03
## Allocationlithium:movement_M10_vmu:phaseweek 2  1.364e-01  1.837e-01  1.196e+03
## Allocationlithium:movement_M10_vmu:phaseweek 3  5.283e-01  1.921e-01  1.198e+03
## Allocationlithium:movement_M10_vmu:phaseweek 4 -5.125e-02  1.954e-01  1.196e+03
##                                                t value Pr(>|t|)    
## (Intercept)                                     -2.088  0.04001 *  
## Allocationlithium                               -2.199  0.02851 *  
## movement_M10_vmu                                 5.472 5.42e-08 ***
## phaseweek 1                                     -3.185  0.00149 ** 
## phaseweek 2                                     -5.689 1.60e-08 ***
## phaseweek 3                                     -5.084 4.29e-07 ***
## phaseweek 4                                     -1.564  0.11798    
## Age                                              0.436  0.66627    
## GenderM                                         -2.068  0.04735 *  
## Allocationlithium:movement_M10_vmu              -2.416  0.01586 *  
## Allocationlithium:phaseweek 1                   -0.067  0.94654    
## Allocationlithium:phaseweek 2                    0.904  0.36624    
## Allocationlithium:phaseweek 3                    3.259  0.00115 ** 
## Allocationlithium:phaseweek 4                    0.187  0.85174    
## movement_M10_vmu:phaseweek 1                    -2.905  0.00374 ** 
## movement_M10_vmu:phaseweek 2                    -5.066 4.71e-07 ***
## movement_M10_vmu:phaseweek 3                    -4.653 3.63e-06 ***
## movement_M10_vmu:phaseweek 4                    -1.494  0.13535    
## Allocationlithium:movement_M10_vmu:phaseweek 1  -0.029  0.97669    
## Allocationlithium:movement_M10_vmu:phaseweek 2   0.743  0.45790    
## Allocationlithium:movement_M10_vmu:phaseweek 3   2.750  0.00605 ** 
## Allocationlithium:movement_M10_vmu:phaseweek 4  -0.262  0.79316    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
pC2 <- plot_model(cB2, type = "pred", terms = c("movement_M10_vmu","Allocation","phase"), show.data = F)+
  scale_color_manual(values=orangeNavyPalette)+
  scale_fill_manual(values=orangeNavyPalette)+
  labs(title = "",
       x = "daytime movement volatility",
       y = "PA volatility")+
  theme(legend.position = "none")
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
pC2$facet$params$nrow=1

cB3 <- lmer(pos_vmu ~ Allocation * time_M10_vmu * phase + Age + Gender + (1|ID),data = df_hlm)
summary(cB3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pos_vmu ~ Allocation * time_M10_vmu * phase + Age + Gender +  
##     (1 | ID)
##    Data: df_hlm
## 
## REML criterion at convergence: 2973.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4155 -0.6696 -0.0367  0.6575  4.1338 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.7574   0.8703  
##  Residual             0.5569   0.7463  
## Number of obs: 1243, groups:  ID, 34
## 
## Fixed effects:
##                                              Estimate Std. Error         df
## (Intercept)                                -2.295e+00  7.187e-01  1.003e+02
## Allocationlithium                           6.977e-01  7.335e-01  5.828e+02
## time_M10_vmu                                4.161e-01  1.244e-01  1.204e+03
## phaseweek 1                                -3.048e-01  6.428e-01  1.196e+03
## phaseweek 2                                -1.291e+00  6.320e-01  1.195e+03
## phaseweek 3                                -1.762e+00  7.224e-01  1.199e+03
## phaseweek 4                                -1.823e+00  7.705e-01  1.200e+03
## Age                                         9.570e-03  1.304e-02  3.023e+01
## GenderM                                    -6.622e-01  3.080e-01  3.004e+01
## Allocationlithium:time_M10_vmu              1.716e-01  1.662e-01  1.202e+03
## Allocationlithium:phaseweek 1              -3.211e+00  9.105e-01  1.199e+03
## Allocationlithium:phaseweek 2              -3.434e+00  8.991e-01  1.199e+03
## Allocationlithium:phaseweek 3              -2.523e+00  9.219e-01  1.199e+03
## Allocationlithium:phaseweek 4              -3.797e+00  9.706e-01  1.200e+03
## time_M10_vmu:phaseweek 1                   -2.202e-03  1.506e-01  1.194e+03
## time_M10_vmu:phaseweek 2                   -1.855e-01  1.449e-01  1.194e+03
## time_M10_vmu:phaseweek 3                   -3.001e-01  1.558e-01  1.196e+03
## time_M10_vmu:phaseweek 4                   -3.314e-01  1.619e-01  1.196e+03
## Allocationlithium:time_M10_vmu:phaseweek 1 -7.456e-01  2.137e-01  1.198e+03
## Allocationlithium:time_M10_vmu:phaseweek 2 -7.923e-01  2.050e-01  1.197e+03
## Allocationlithium:time_M10_vmu:phaseweek 3 -6.578e-01  2.042e-01  1.196e+03
## Allocationlithium:time_M10_vmu:phaseweek 4 -9.397e-01  2.101e-01  1.197e+03
##                                            t value Pr(>|t|)    
## (Intercept)                                 -3.193 0.001881 ** 
## Allocationlithium                            0.951 0.341889    
## time_M10_vmu                                 3.344 0.000852 ***
## phaseweek 1                                 -0.474 0.635452    
## phaseweek 2                                 -2.042 0.041344 *  
## phaseweek 3                                 -2.438 0.014893 *  
## phaseweek 4                                 -2.366 0.018156 *  
## Age                                          0.734 0.468848    
## GenderM                                     -2.150 0.039726 *  
## Allocationlithium:time_M10_vmu               1.032 0.302095    
## Allocationlithium:phaseweek 1               -3.526 0.000437 ***
## Allocationlithium:phaseweek 2               -3.819 0.000141 ***
## Allocationlithium:phaseweek 3               -2.737 0.006287 ** 
## Allocationlithium:phaseweek 4               -3.912 9.67e-05 ***
## time_M10_vmu:phaseweek 1                    -0.015 0.988340    
## time_M10_vmu:phaseweek 2                    -1.281 0.200573    
## time_M10_vmu:phaseweek 3                    -1.926 0.054343 .  
## time_M10_vmu:phaseweek 4                    -2.047 0.040869 *  
## Allocationlithium:time_M10_vmu:phaseweek 1  -3.490 0.000501 ***
## Allocationlithium:time_M10_vmu:phaseweek 2  -3.865 0.000117 ***
## Allocationlithium:time_M10_vmu:phaseweek 3  -3.221 0.001311 ** 
## Allocationlithium:time_M10_vmu:phaseweek 4  -4.474 8.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 22 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
pC3 <- plot_model(cB3, type = "pred", terms = c("time_M10_vmu","Allocation","phase"), show.data = F)+
  scale_color_manual(values=orangeNavyPalette)+
  scale_fill_manual(values=orangeNavyPalette)+
  labs(title = "",
       x = "daytime movement onset time volatility",
       y = "PA volatility")+
  theme(legend.position = "none")
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
pC3$facet$params$nrow=1
pC1+pC2+pC3+
  plot_layout(nrow = 3)+
  plot_annotation(tag_levels = c("A"))

ggsave("../figs/FIG5_change_slope.png",width = 7 ,height=8,dpi = 300)